{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1：$\\theta$ 是三⻆视差，其测量服从 $N(1,0.09)$，⽤随机抽样方法给出距离 $d=1/\\theta$ 的误差分布\n",
    "\n",
    "2：求 $\\theta$ 的平均值（数学期望），⽅差；$d$ 的数学期望和⽅差"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from scipy import integrate\n",
    "\n",
    "matplotlib.rc('xtick', labelsize=12) \n",
    "matplotlib.rc('ytick', labelsize=12) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 5)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAF+CAYAAADQoS0WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYXVWd7vHvLwNJyAQkYSYMkgQBO1GKBA0YBBRBaaDRCwoIdisKF23FuS8ojba2XkCviiDdIEMrDsgoCIqITIoUaoAICYMhkIFUBEICmbPuH6sKKkUNp5Kqs8/Z5/t5nvPsql37VL1YkrystfdakVJCkiRJtW1A0QEkSZLUM0ubJElSHbC0SZIk1QFLmyRJUh2wtEmSJNUBS5skSVIdsLRJkiTVAUubJElSHbC0SZIk1QFLmyRJUh0YVHSAvjZ27Ni0yy67FB1DkiSpRw888MCSlNK4Sq6temmLiNOBk4E3AFellE7u5trdgG8DM4BVwKUppc929/132WUXmpub+yyvJElSf4mIpyq9tojp0QXAV4BLu7soIjYDfg3cDmwL7Aj8T7+nkyRJqkFVH2lLKV0DEBFN5CLWlZOBBSml89ude7Afo0mSJNWsWn4QYT9gbkT8MiKWRMQdEfGGokNJkiQVoZZL247AceR72rYHbgKub5023UBEnBIRzRHR3NLSUuWYkiRJ/a+WS9sK4O6U0i9TSquBc4ExwOs7XphSujil1JRSaho3rqIHMCRJkupKLZe2B4FUdAhJkqRaUPXSFhGDImIoMBAYGBFDI6KzByL+B9gvIg6JiIHAJ4AlwCNVjCtJklQTihhpO5M89fl54ITWj8+MiPERsTwixgOklGa3fv0i4HngSOAfW6dKJUmSGkqkVK4ZyKampuTiupIkqR5ExAMppaZKrq3le9okSZLUytImSZJUByxt6jMvvAC33gr33190EkmSysfSpj5x880waRK8850wdSocdhg8/3zRqSRJKg9LmzbZQw/Be94D220Ht9wC550Ht98ORx4JK1cWnU6SpHKo+obxKpf16+GEE2CLLXJh23ZbOPRQ2HFHOPZYOOcc+OpXi04pSVL9c6RNm+Tqq+HBB/Po2rbbvnr+f/0vOOkk+L//F2bNKi6fJEllYWnTRkspj6TtuWcuaR2dey5svjmcfXbVo0mSVDqWNm20P/4xj6KdcQYMHPjar48dCx/7GPz85/DXv1Y/nyRJZWJp00a7/HIYNgze+96ur/nEJ/I1559fvVySJJWRpU0bZfVq+PGP4eijYdSorq8bOxbe/3646ipYurR6+SRJKhtLmzbKPffkddi6G2Vr85GPwMsvw5VX9n8uSZLKytKmjXLzzTB4MBxySM/XNjXBlCl5OlWSJG0cS5s2ys03w4wZMGJEZdefcAI0N8Njj/VvLkmSysrSpl6bNy8/DXr44ZW/59hjISLf2yZJknrP0qZeu/POfDzooMrfs+OOcMAB8LOf9U8mSZLKztKmXrvrLhg9Gvbeu3fvO/JIePhhmDu3X2JJklRqljb12l13wfTpnS+o250jjsjHG2/s+0ySJJWdpU29smQJPPII7L9/7987YQJMmgS/+EXf55IkqewsbeqVP/whH6dP37j3H3EE3HEHLFvWZ5EkSWoIljb1SnNzfgr0TW/auPe/+915N4Vf/apvc0mSVHaWNvXKAw/AHntUvj5bR9Onw5Zbel+bJEm9ZWlTrzzwAOyzz8a/f9AgOOwwuOkmWL++73JJklR2ljZVbOHC/NqU0gbwznfmBxoefLBvckmS1AgsbarYAw/k46aWtre9LR9vv33Tvo8kSY3E0qaKtY2MTZ68ad9nxx1h4kRLmyRJvWFpU8Uefhh23hlGjdr073XQQfC738GaNZv+vSRJagSWNlXs4Yd7v3VVVw4+GJYvz0uISJKknlnaVJE1a+DRR/uutB14YD46RSpJUmWqXtoi4vSIaI6IVRFxWYXvuT0iUkQM6ud46sJjj+Xi1lelbezYfG+cpU2SpMoUMdK2APgKcGklF0fE8YBlrWAPP5yPfVXaIN/Xds89sHJl331PSZLKquqlLaV0TUrpOuDvPV0bEaOBLwGf7fdg6tasWTBgQN4Noa8cdBCsWvXqfqaSJKlrtX5P21eBC4FF3V0UEae0Trk2t7S0VCdZg5k9G3bZBYYO7bvv2bbp/N139933lCSprGq2tEVEEzAd+E5P16aULk4pNaWUmsaNG9f/4RrQ7NkwaVLffs8tt8zTrZY2SZJ6VpOlLSIGAN8D/jWltLboPI0uJZgzJy+I29cOOADuvRfWrev77y1JUpnUZGkDRgFNwE8iYhFwf+v5ZyLigOJiNab58+Hll/t+pA1g//1h2TJ46KG+/96SJJVJ1Z/KbF22YxAwEBgYEUOBtR1G1JYC27f7fCfgj8A+gDetVdns2fnYX6UN8hTplCl9//0lSSqLIkbazgRWAJ8HTmj9+MyIGB8RyyNifMoWtb14tag9m1JaXUDmhjZnTj72x/To+PGw007e1yZJUk+qPtKWUjobOLuLL4/o4j1zgeifROrJ7NkwfDjssEP/fP/998/7kKYE4W9ZkqRO1eo9baohjz8Ou+/ef4Vq//1hwQJ46qn++f6SJJWBpU09evJJ2G23/vv+bfe13XVX//0MSZLqnaVN3Vq/Ppe2172u/37GXnvByJHujCBJUncsberWggV5q6n+LG0DB8LUqZY2SZK6Y2lTt558Mh/7c3oUYNo0ePDBvB6cJEl6LUubuvXEE/nYnyNtAPvtB2vXwp/+1L8/R5KkemVpU7eeeCJPX44f378/Z9q0fHSKVJKkzlna1K0nnsiFbfDg/v05W28Nu+4K993Xvz9HkqR6ZWlTt+bOzWWqGvbbz5E2SZK6YmlTt+bNg513rs7PmjYNnnkmb1AvSZI2ZGlTl1avhoUL+/9+tjb77ZePTpFKkvRaljZ1af78vB9otUrblCmw2WZOkUqS1BlLm7r09NP5WK3SNmQIvPGNljZJkjpjaVOX5s3Lx2qVNshTpM3Nec02SZL0KkubutRW2nbaqXo/c9o0WLECHnqoej9TkqR6YGlTl+bNg3HjYNiw6v3MtocRnCKVJGlDljZ1ad686k6NAuyyS15o1ydIJUnakKVNXSqitEXkKVJH2iRJ2pClTZ1KCZ56qvqlDWDqVJg9G5Yurf7PliSpVlna1KmlS2H58uJKG+SnSCVJUmZpU6eKWO6jTVNTPv7xj9X/2ZIk1SpLmzpVZGnbaiuYMMHSJklSe5Y2darI0gZ5itTSJknSqyxt6tS8eXkf0K23LubnT50KCxbk/U8lSZKlTV2YNy/vhDCgoP+HtD2M4GibJEmZpU2dKmKNtvamTIFBgyxtkiS1sbSpU0WXtqFDYfJkS5skSW0sbXqNdevy/WQ77lhsjqlT81pt69cXm0OSpFpQ9dIWEadHRHNErIqIy7q57qSIeCAiXoyIZyLiGxExqIpRG1ZLSy5u229fbI6pU+HFF2HOnGJzSJJUC4oYaVsAfAW4tIfrNgc+AYwFpgEHA5/u32iCPMoGxZe2fffNR6dIJUkqoLSllK5JKV0H/L2H6y5MKd2VUlqdUpoP/BCYXpWQDW7hwnwsurTtsQeMGGFpkyQJ6uuetrcCs4oO0QhqZaRt4MC8pZWlTZKkOiltEfFBoAk4t4uvn9J6n1xzS0tLdcOV0IIFEAHbbFN0knxf21/+AqtWFZ1EkqRi1Xxpi4ijgP8EDkspLensmpTSxSmlppRS07hx46obsIQWLIBx42Dw4KKT5NK2Zg3MnFl0EkmSilXTpS0i3gn8F3BESumhovM0igULip8abePOCJIkZVVfQqN12Y5BwEBgYEQMBdamlNZ2uO4g8sMHR6eU/Cu7imqptO24I2y7raVNkqQiRtrOBFYAnwdOaP34zIgYHxHLI6JtHf6zgNHAza3nl0fELwvI23AWLqyd0haRR9ssbZKkRlf1kbaU0tnA2V18eUS7695WjTza0Nq18OyztVPaIJe2G26AF16ALbYoOo0kScWo6XvaVH2LF+dto7bbrugkr2q7r625udgckiQVydKmDdTKGm3tNTXlo1OkkqRGZmnTBmqxtG25JUycCPffX3QSSZKKY2nTBmqxtIEPI0iSZGnTBhYuhAEDYOuti06yoX33zYVy/vyik0iSVAxLmzawYEEubIOq/lxx91xkV5LU6Cxt2kAtLazb3pQpuUha2iRJjcrSpg3UamkbOhQmT7a0SZIal6VNG6jV0gZ5ivT++/M6cpIkNRpLm16xZk1eXLeWS9uyZTB7dtFJJEmqPkubXvHss/lYy6UNnCKVJDUmS5te0bZG27bbFpujK5MmwciRljZJUmOytOkVbSNttVraBg7MW1pZ2iRJjcjSple0lbZttik2R3emToWZM2HVqqKTSJJUXZY2vaKttNXabgjtTZ2aH5iYObPoJJIkVZelTa9YtAhGj85rotUqH0aQJDUqS5te8eyztXs/W5sddsgZLW2SpEZjadMrnn22tu9nA4jIo22WNklSo7G06RX1UNogl7bZs+GFF4pOIklS9Vja9Ip6Km0Azc3F5pAkqZosbQLyEhovvFAfpa2pKR+dIpUkNRJLm4DaX1i3vS23hIkTLW2SpMZiaRNQHwvrtjd1Ktx3H6RUdBJJkqrD0iagPkvbokUwf37RSSRJqg5Lm4D6LG3gFKkkqXFY2gTUX2mbPBkGD4b77y86iSRJ1WFpE5BL26hRtb2FVXtDh+bi5kibJKlRWNoE5PvD6uHJ0famTs0jbevXF51EkqT+V/XSFhGnR0RzRKyKiMt6uPaTEbEoIpZGxKURMaRKMRtOvSys297UqbBsWd4dQZKksitipG0B8BXg0u4uiohDgc8DBwO7ALsB/97f4RpVPZa2fffNR6dIJUmNoOqlLaV0TUrpOuDvPVx6EnBJSmlWSul54MvAyf2dr1HVY2mbNAlGjrS0SZIaQy3f07YXMLPd5zOBbSJiTEF5Smv1anj++forbQMH5i2tLG2SpEZQy6VtBLC03edtH4/seGFEnNJ6n1xzS0tLVcKVyeLF+VhvDyJAvq9t5kxYubLoJJIk9a9aLm3LgVHtPm/7eFnHC1NKF6eUmlJKTePGjatKuDJZtCgf622kDXJpW7MmFzdJksqslkvbLGByu88nA8+mlHq6F069VG8L67bnzgiSpEZRxJIfgyJiKDAQGBgRQyNiUCeXXgH8S0TsGRFbAmcCl1UxasOo59K2ww6w3XaWNklS+RUx0nYmsIK8nMcJrR+fGRHjI2J5RIwHSCndAnwD+C3wVOvrSwXkLb22e9q23rrYHBsjIo+2WdokSWXX2QhXv0opnQ2c3cWXR3S49nzg/H6O1PBaWmD4cNh886KTbJypU+H66+GFF2CLLYpOI0lS/6jle9pUJYsXQz0/v9F2X1tzc7E5JEnqT5Y20dJSn1OjbZqa8tEpUklSmVnaVPcjbVtskXdHsLRJksrM0qa6H2mDvA/pffdBSkUnkSSpf1jaGlxKubTV80gb5PvaFi2C+fOLTiJJUv+wtDW4Zctg1ar6H2lzkV1JUtlZ2hpc21at9T7SNnkyDB5saZMklZelrcG1Laxb76Vt6NBc3CxtkqSysrQ1uLaRtnqfHoU8RdrcDOvWFZ1EkqS+Z2lrcGUZaYNc2pYtg9mzi04iSVLfs7Q1uLLc0wYwbVo+/uEPxeaQJKk/WNoaXEsLjBgBw4YVnWTTTZwIW20F99xTdBJJkvqepa3B1ftuCO0NGADTp8PddxedRJKkvmdpa3Bl2A2hvenTYc6cV6d9JUkqC0tbgyvTSBvk0gZw773F5pAkqa9Z2hpc2Ubamppgs82cIpUklY+lrYGVZd/R9oYOzcXNhxEkSWVjaWtgL74Iq1eXa6QN8hRpczOsWFF0EkmS+o6lrYGVaY229qZPhzVrcnGTJKksLG0NrEy7IbT3lrfko1OkkqQysbQ1sDLtO9reuHEwaZKlTZJULpa2BlbWkTbIU6T33APr1xedRJKkvmFpa2BlvacNYP/94fnn4dFHi04iSVLfsLQ1sJYWGDkyL5NRNm2L7DpFKkkqC0tbA1u8uHz3s7WZMCGPIN51V9FJJEnqG5a2Bla2hXXbi4ADDoA77yw6iSRJfcPS1sDKtu9oRzNmwFNP5ZckSfXO0tbAyrbvaEcHHpiPv/tdoTEkSeoTFZe2iNgmIk6MiHMj4r9bjydGxLa9+YERsVVEXBsRL0XEUxHx/i6uGxIRF0XEsxHxXETcGBE79OZnqWtl3He0o733hq22gjvuKDqJJEmbrsfSFhGvj4irgb8CJwKDgUWtxxOBWRFxdUTsWeHPvABYDWwDHA9cGBF7dXLdvwJvBv4B2B54AfhOhT9DPVi6NG/1VOaRtgED8n1tjrRJkspgUAXXXAacCxyfUlrV8YsRsRlwJHAJuWR1KSKGA8cAe6eUlgN3R8QN5PL3+Q6X7wrcmlJ6tvW9PwbOryCvKlDmNdramzEDrr8ennkGdtyx6DSSJG28HkfaUkrTUko/66ywtX59devXuy1srSYC61JKc9qdmwl0NtJ2CTA9IraPiM3Jo3K/rOBnqAJl3g2hPe9rkySVRSXTo9+OiA9HxLTWkbJNMQJY2uHcUmBkJ9fOAeYB84EXgdcD53SR8ZSIaI6I5pa2ISR1q6z7jnb0D/8Ao0d7X5skqf5V8iBCC3AQ8F/A4oh4PCKuiYh/j4hjImJiRESFP285MKrDuVHAsk6uvRAYCowBhgPX0MVIW0rp4pRSU0qpaVzZh476SKOMtA0c6H1tkqRyqGR69MsppfellP4B+GfyQwj3A6OB84BHyGWsEnOAQRExod25ycCsTq6dDFyWUnqudWr2O8DUiBhb4c9SNxrlnjbIU6SPPQYLFxadRJKkjdfbddrOAw5LKX0tpfQJYALwA+Bblbw5pfQSecTsnIgYHhHTyQ8xXNnJ5fcDH4iI0RExGDgNWJBSWtLLzOpESwuMGgVDhhSdpP/NmJGPjrZJkupZb0vbOvIIGwAppTXAGeQnQit1GjAMWAxcBZyaUpoVEQdERPsRu08DK4HHyFO0hwNH9zKvulDmfUc7mjIFRo70vjZJUn2rZMmP9i4CroqIE1NKc1vPbQNsV+k3SCk9BxzVyfm7yA8qtH3+d/ITo+oHZV9Yt71Bg7yvTZJU/3o70vafwG+AhyLivoi4Efgj8KM+T6Z+VfZ9Rzs68EB49FHva5Mk1a9elbaUnU1e+PZ84BbgvSmlU/shm/pR2fcd7eiQQ/LxN78pNockSRurknXaPh4RG9yunlJaklL6SUrpgpTSba37hH68/2KqLzXCvqMdTZ4MY8bAbbcVnUSSpI1TyT1t2wKPR8TNwO+A2eR11UaSdziYQX5I4Ir+Cqm+9cILsHZtY420DRgABx+cS1tKUPHKgpIk1YhKpkfHAW8kP8X5L+QFbh8Gbiav2zYbeGNK6cz+Cqm+1UhrtLV3yCEwfz7Mnl10EkmSeq+SkbZ/Sil9GDi39anRBhqfKae23RAaaaQNXr2v7bbbYI89is0iSVJvVTLS9kBEfDciZpAfQFCda9SRtl13hd128742SVJ9qqS0vR9YBXwbGB4RCyPiloj4ekS8PyL2ioiB/RtTfalR9h3tzCGHwG9/m+/pkySpnlSy9+iSlNKnUkqTyQ8gHARc3vrlDwC3Ufneo6oBjTrSBrm0vfgiNDcXnUSSpN7p7Y4I41q3rnqEvAUVABGxTZ+mUr9qaYHRo2GzzYpOUn1ve1t+cvS222C//YpOI0lS5Xq7uO6aLs4/2zdxVA2NtO9oR2PHwhvf6H1tkqT609ttrFQCjbawbkeHHAL33gsvvVR0EkmSKmdpa0CNtu9oR4ccAmvWwJ13Fp1EkqTKWdoaUKOPtO2/PwwdCrfcUnQSSZIqZ2lrMOvX59K2TQM/OjJsGBx0ENx8c9FJJEmqnKWtwTz/PKxb17gPIrQ5/HB4/HF47LGik0iSVBlLW4Np5DXa2jvssHz85S+LzSFJUqUsbQ2mUfcd7Wi33WDSJKdIJUn1w9LWYCxtrzr8cLjjDnj55aKTSJLUM0tbg7G0veqww2DVqrwXqSRJtc7S1mAWL87bOI0ZU3SS4r31rbD55k6RSpLqg6WtwbS0wFZbwaDe7jpbQkOG5IV2b74ZUio6jSRJ3bO0NZhG3ne0M4cdBnPnwuzZRSeRJKl7lrYGY2nbkEt/SJLqhaWtwVjaNrTzzrDXXnDjjUUnkSSpe5a2BmNpe60jj8ybxz/3XNFJJEnqmqWtgaxdm4tJo++G0NFRR+WtvW66qegkkiR1zdLWQJYsyUdH2ja0zz6www5w3XVFJ5EkqWtVL20RsVVEXBsRL0XEUxHx/m6ufVNE3BkRyyPi2Yj412pmLRsX1u3cgAF5ivSWW2DFiqLTSJLUuSJG2i4AVgPbAMcDF0bEXh0vioixwC3A94ExwO7Ar6qYs3QsbV076qi8ndVttxWdRJKkzlW1tEXEcOAY4KyU0vKU0t3ADcCJnVx+BnBrSumHKaVVKaVlKaVHqpm3bFpa8tF72l5rxgwYPdopUklS7ar2SNtEYF1KaU67czOB14y0AfsBz0XEvRGxOCJujIjxVUlZUo60dW2zzeBd74IbbsgPJUiSVGuqXdpGAEs7nFsKjOzk2h2Bk4B/BcYDfwOu6uybRsQpEdEcEc0tbcNJeo3Fi/P2VVtsUXSS2nTUUflhjXvvLTqJJEmvVe3SthwY1eHcKGBZJ9euAK5NKd2fUloJ/DvwlogY3fHClNLFKaWmlFLTOOf+urR4cZ4aHeAzw5165zvziJtTpJKkWlTtv77nAIMiYkK7c5OBWZ1c+yDQfhvvto+jn7KVngvrdm/kyLyB/LXXuoG8JKn2VLW0pZReAq4BzomI4RExHTgSuLKTy38AHB0RUyJiMHAWcHdK6YXqJS6XlhYfQujJ0UfD3/4Gf/5z0UkkSdpQERNlpwHDgMXke9ROTSnNiogDImJ520UppduBfwNuar12d6DLNd3UM0faenb00fm+v5/8pOgkkiRtqOqlLaX0XErpqJTS8JTS+JTSj1rP35VSGtHh2gtTSjuklLZMKR2RUnq62nnLxNLWszFj4O1vz6XNKVJJUi3xlvQGsWIFLFtmaavEccfBU0/BffcVnUSSpFdZ2hpE20oolraeHXlkforUKVJJUi2xtDUId0Oo3OjRcPjh8NOfwvr1RaeRJCmztDUId0PonWOPhQUL4O67i04iSVJmaWsQlrbeefe7Ydgwp0glSbXD0tYgLG29M2IEHHEE/OxnsHZt0WkkSbK0NYzFi2HoUBg+vOgk9ePYY/O9gL/9bdFJJEmytDWMlpY8yhZuAlaxww/PDyVc2dl+HZIkVZmlrUG4sG7vDR2a12z7+c/zGneSJBXJ0tYgLG0b5+ST4eWX4eqri04iSWp0lrYGYWnbONOmwcSJcNllRSeRJDU6S1sDSCnf0+bCur0XkUfb7rwTnnyy6DSSpEZmaWsAy5fDypWOtG2sE0/M5e2KK4pOIklqZJa2BuAabZtmxx3hkEPg8svd1kqSVBxLWwNoK21Oj268k0+GuXPzNKkkSUWwtDWARYvycdtti81Rz446CkaNgh/8oOgkkqRGZWlrAJa2Tbf55vC+98FPfwrPPVd0GklSI7K0NYBFi/KN9E6PbppTT80PdFx+edFJJEmNyNLWABYtyoVt0KCik9S3yZPhzW+Giy7Ky6hIklRNlrYGsGgRbLdd0SnK4dRTYc4cuP32opNIkhqNpa0BLFrk/Wx95b3vhTFj4MILi04iSWo0lrYGYGnrO0OHwgc/CNddBwsWFJ1GktRILG0ll5Klra995COwbh38938XnUSS1EgsbSX3/POwerWlrS/tvju84x1w8cWwdm3RaSRJjcLSVnKu0dY/Tj8d5s+Hq68uOokkqVFY2krO0tY/3vUumDgRzj3X5T8kSdVhaSs5S1v/GDAAzjgDHnjA/UglSdVhaSs5S1v/+cAHYOzYPNomSVJ/s7SV3KJFMGQIjB5ddJLyGTYM/vf/hl/8Ah59tOg0kqSyq3ppi4itIuLaiHgpIp6KiPf3cP1mEfFoRDxTrYxl0rbcR0TRScrptNPy2m3nn190EklS2RUx0nYBsBrYBjgeuDAi9urm+s8Ai6sRrIxco61/bb01nHQSXHEFPPts0WkkSWVW1dIWEcOBY4CzUkrLU0p3AzcAJ3Zx/a7ACcDXqpeyXCxt/e+Tn8xr4f2//1d0EklSmVV7pG0isC6lNKfduZlAVyNt3wH+DVjR38HKauFCS1t/mzQp70n6ne/A3/9edBpJUllVu7SNAJZ2OLcUGNnxwog4GhiUUrq2p28aEadERHNENLe0tPRN0hJYswaWLLG0VcNZZ8Hy5fDNbxadRJJUVtUubcuBUR3OjQKWtT/ROo36DeBjlXzTlNLFKaWmlFLTuHHj+iRoGSxuvRPQ0tb/9t47j7Z9+9vw3HNFp5EklVG1S9scYFBETGh3bjIwq8N1E4BdgLsiYhFwDbBdRCyKiF2qkLMUXKOtus46C5Ytc7RNktQ/qlraUkovkQvYORExPCKmA0cCV3a49GFgJ2BK6+tDwLOtHz9dvcT1zdJWXW94A7znPY62SZL6RxFLfpwGDCMv43EVcGpKaVZEHBARywFSSmtTSovaXsBzwPrWz9cVkLkuWdqq74tfhBdfhG99q+gkkqSyqXppSyk9l1I6KqU0PKU0PqX0o9bzd6WURnTxnjtSSjtWN2n9W7gwH7fZptgcjaRttO2b33TdNklS33IbqxJ75hkYMyZvt6Tq+Y//gJUr4Zxzik4iSSoTS1uJzZ8POzo+WXUTJ8JHPgLf/z7Mnl10GklSWVjaSuyZZyxtRfniF/MI5xe+UHQSSVJZWNpKzNJWnK23hs99Dq69Fu65p+g0kqQysLSV1MqVeTeEHXYoOknjOuMM2H57+MxnIKWi00iS6p2lraQWLMhHR9qKs/nm8OUvw+9/Dz/6UdFpJEn1ztJWUs88k4+WtmKdfDJMnQqf+hQs7bjrriRJvWBpKylLW20YMAC+9728D+wXv1h0GklSPbO0lVRbafOetuLMln8HAAASh0lEQVTtsw989KPw3e/CX/5SdBpJUr2ytJXU/PkwciSMGlV0EkFecHfMGDjtNFi/vug0kqR6ZGkrKZf7qC1bbgnf+EZ+KOHSS4tOI0mqR5a2krK01Z4PfABmzMgPJTz9dNFpJEn1xtJWUs884/1stWbAgDzKtnYtfPjDrt0mSeodS1sJrV0LixY50laLdtsNvv51uPVW+MEPik4jSaonlrYSWrQo3+xuaatNp50GBx4In/yk06SSpMpZ2krINdpq24ABcMklsG4d/Mu/+DSpJKkylrYSco222rfbbnDeefDrX+ejJEk9sbSV0Pz5+ehIW2075RR4z3vg3/4N/vCHotNIkmqdpa2EnnkGhgzJi7mqdkXAf/1XLtfvex+88ELRiSRJtczSVkJta7RFFJ1EPdliC7jqqvw7+9CHXAZEktQ1S1sJzZvn1Gg92W+/vM3Vz38O3/xm0WkkSbXK0lZCc+fCrrsWnUK98elPwzHHwGc+A7/6VdFpJEm1yNJWMqtWwYIFlrZ6M2AAXHYZ7LUXHHssPP540YkkSbXG0lYy8+bl4y67FBpDG2HECLj++lzg/vEf4cUXi04kSaollraSmTs3Hy1t9WnXXeFnP4M5c/KI25o1RSeSJNUKS1vJWNrq30EHwUUXwS23+ESpJOlVg4oOoL41dy4MGgTbb190Em2KD30o35v4pS/lnS2++tWiE0mSimZpK5m5c2GnnXJxU30766xc3L72NdhuO/jYx4pOJEkqUtWnRyNiq4i4NiJeioinIuL9XVz3mYh4OCKWRcTfIuIz1c5aj/72N9h556JTqC9EwAUXwFFHwcc/nndPkCQ1riLuabsAWA1sAxwPXBgRe3VyXQAfALYE3gmcHhHHVS1lnXriCXjd64pOob4ycCD8+Mdw2GF5r9JLLy06kSSpKFUtbRExHDgGOCultDyldDdwA3Bix2tTSt9IKf0ppbQ2pTQbuB6YXs289ebFF2HxYth996KTqC8NGQLXXAPveEe+1+3yy4tOJEkqQrVH2iYC61JKc9qdmwl0NtL2iogI4ABgVj9mq3tPPJGPEyYUm0N9b+hQuO46OPhg+OAH4fvfLzqRJKnaql3aRgBLO5xbCozs4X1nk7P+oLMvRsQpEdEcEc0tLS2bHLJeta2i70hbOQ0blhffPeww+OhH836lLgciSY2j2qVtOTCqw7lRwLKu3hARp5PvbXtXSmlVZ9eklC5OKTWllJrGjRvXZ2HrTVtp85628tp88zzidvzxcOaZcMYZsH590akkSdVQ7YUh5gCDImJCSumx1nOT6WLaMyL+Gfg88NaU0jNVyli3Hn8ctt02b4ek8ho8GK64AsaOhW99CxYuhB/8II/ESZLKq6ojbSmll4BrgHMiYnhETAeOBK7seG1EHA98FXh7SunJauasV48/7tRooxgwAL75Tfj61+GnP4UZM/KabpKk8ipiyY/TgGHAYuAq4NSU0qyIOCAilre77ivAGOD+iFje+rqogLx147HHLG2NJAI++1m49lr4619h6lR44IGiU0mS+kvVS1tK6bmU0lEppeEppfEppR+1nr8rpTSi3XW7ppQGp5RGtHt9tNp568XSpXmabI89ik6iajvySLjnnrym2/77wyWX+ICCJJWRG8aXxCOP5OPrX19sDhVj8mS4//5c2j70ITj5ZHjppaJTSZL6kqWtJNpK2557FptDxdl6a7jlFjj7bLjySpg2DR56qOhUkqS+Ymkrib/+Na+cv+uuRSdRkQYOhC99CW69FZYsgaYm+MY3YN26opNJkjaVpa0kHnkEJk3Kf2lLb397HmV797vhc5+DAw98dccMSVJ9srSVxF//6tSoNjRuHFx9dZ4qfegh2Htv+OpXYfXqopNJkjaGpa0EXn4Z5s71IQS9VgSccALMmgVHHAH/5//AlCnwu98VnUyS1FuWthJ48MG8xMPkyUUnUa3aYYe8CO9NN8HKlXm69KSTYP78opNJkiplaSuBv/wlH6dMKTaHat/hh8PDD8MXvgA//jFMmJBH35YuLTqZJKknlrYS+POfYcstYfz4opOoHmy+eb637dFH4eij88e77w7f+Y73u0lSLbO0lcCf/5xH2SKKTqJ6suuu8MMfQnMzvOEN8PGP5/L23e/CihVFp5MkdWRpq3Nr1+YnA9/4xqKTqF7tsw/85jd5Yd7x4+FjH8uF7txzYfnynt8vSaoOS1ude+SRfGO5pU2bIgIOPRTuugvuuCOPvH3mM7nEffaz8NRTRSeUJFna6tzvf5+P++1XbA6VQwTMmAG//nX+/9bBB8P558Nuu8E//VMudG5GL0nFsLTVuXvvzYuovu51RSdR2ey3H/zsZ/Dkk3m07c474W1vy+sBfv3rsHBh0QklqbFY2urcvffCW97iQwjqP+PHw9e+Bk8/DZdemv8j4fOfh512yttkXX11XuBZktS/LG11rKUFHnsslzapvw0bBh/8YL7vbfbsPPr25z/De98LW28Nxx0H11zjk6eS1F8sbXXsrrvy0dKmaps4Ma/vNm9efvL0hBPg9tvhmGPySNz73gc/+Qk8/3zRSSWpPCxtdexXv4KRI2HatKKTqFENHAgHHQQXXQQLFsBtt8Hxx+fjccfB2LFwwAF5erVtuzVJ0saJVLI/RZuamlJzc3PRMfpdSvmJvsmT4brrik4jbWjdOvjjH+Hmm/PrT3/K57ffPj/M8La35f1Pd9vN+zElNbaIeCCl1FTRtZa2+jRnDkyaBN/7Hpx6atFppO4tWJAX77311rxsyOLF+fxOO+XydsAB+WnVPffMo3eS1Ch6U9oG9XcY9Y9f/CIfDz202BxSJbbfHv75n/Mrpbwo9B135Nctt8CVV+brRoyAfffNBW7atPzadtsik0tS7XCkrU7tu2/+y68B/lFVcinB44/DfffBH/6QXzNn5i3aALbbLu+tO3lyfk2ZAhMmOCInqRwcaSu5xx7LZe3cc4tOIm26iFzCJkzIT6FCXjbkT3/KRe4vf8kl7te/frXIDRsGe+2VF/rdY49XX697HQwZUtw/iyT1J0tbHWqbSjr22GJzSP1l2DCYPj2/2qxenadVZ87MRe6hh/L0atu/D5BH33bdNRe4iRPzx22vXXaBzTev9j+JJPUdp0frzIoVeYX6/faDG28sOo1UvOXL84M5jz664euxx2Dlyg2v3XbbDYvczjvDDju8+hozxqdZJVWX06MldsUVsGQJfOpTRSeRasOIEfCmN+VXeynBs8/C3/624evJJ/P2bz/5SV6apL0hQ/JDE9tvv2GZ2267vOvDuHH5OHYsDB5cvX9GSQJH2urKiy/mZT523hl+/3tHBKRNsXZt3vR+/vzuX13tq7rllq+WuHHjNix0W2654WuLLfJx2DD/vZW0IUfaSuqss/LIwfXX+we/tKkGDcrrxO20U9fXpARLl8KiRXmv38WLOz/OmQP33JNHwdev7/r7bbbZa4tc22vkyA1fI0Z0/bkPW0iNqeqlLSK2Ai4B3gEsAb6QUvpRJ9cF8J/Ah1pPXQJ8LpVtaLBCP/whfPvbcPrpMHVq0WmkxhCRy9UWW+SHG3qybh288ELec7Xt1fHz9q/Fi2H27Pzxiy++drq2K4MHb1jiRozIo3ibb975sdJzQ4fmQtj22myzXG79j0SpNhQx0nYBsBrYBpgC3BQRM1NKszpcdwpwFDAZSMCvgSeBi6qYtXDr18MFF8AnPwkzZsB55xWdSFJXBg7MDzOMGdP796YEq1bBsmX54Yply159dfd528crVuSR+BUr8pRu+2PbUikbI+LVAte+zPX248GDcwFsf+zsXF8cBw3Kv4v2rwEDLJ+qf1UtbRExHDgG2DultBy4OyJuAE4EPt/h8pOA81JKz7S+9zzgwzRAaVu7FubNg9tvzxtxP/AAvPvdcNVV+Q8/SeUTkUe6hg7N98f1pTVrcnnrrNB1/Hj16lwe245tr/afd/bxSy/lEcOurlmzJv/Z1t30cX+LeG2Zayt0lZzb1PMR+fO2Atndx5v69f66NmLDV9v/rpWcr+a1tZir7fP2x0G9bGHVHmmbCKxLKc1pd24mMKOTa/dq/Vr76/bq6Qc8+ii8+c35v1rhtcfOzlXytWq9/+WX8/0zbX+wTZwIl10GJ56Y/4WRpN5qG9UaNaroJPnPtrVrXy1x/XVct27D1/r1rz3X2/PdXbtmTc/Xp5Rf69fnV08fV3ptY940VA6TJvXu+mqXthHA0g7nlgIjK7h2KTAiIqLjfW0RcQp5OpWhQ9/wyh9MXTXbjf3apr6/ku89ZEheYmD8+Lzv4hve4JC+pPIYMCDPGDhr0Lf6qgBW+vX2gw4dX52dr+a1tZir7fOOxy22gJNOqvz3XO3Sthzo+N96o4BlFVw7Clje2YMIKaWLgYshL/lx6619E1aSpHrQfupX9aU3pa3aE25zgEERMaHduclAx4cQaD03uYLrJEmSSq+qpS2l9BJwDXBORAyPiOnAkcCVnVx+BXBGROwQEdsDnwIuq1pYSZKkGlLEre2nAcOAxcBVwKkppVkRcUBELG933feBG4GHgIeBm1rPSZIkNZyqr9OWUnqOvP5ax/N3kR8+aPs8AZ9tfUmSJDU0F5GQJEmqA5Y2SZKkOmBpkyRJqgOWNkmSpDpgaZMkSaoDljZJkqQ6YGmTJEmqA5Y2SZKkOmBpkyRJqgORNx4oj4hYBswuOoc22lhgSdEhtFH83dU3f3/1zd9f/ZqUUhpZyYVV38aqCmanlJqKDqGNExHN/v7qk7+7+ubvr775+6tfEdFc6bVOj0qSJNUBS5skSVIdKGNpu7joANok/v7ql7+7+ubvr775+6tfFf/uSvcggiRJUhmVcaRNkiSpdCxtkiRJdaA0pS0itoqIayPipYh4KiLeX3QmVSYiTo+I5ohYFRGXFZ1HlYuIIRFxSeu/c8si4s8RcVjRuVS5iPifiFgYES9GxJyI+FDRmdQ7ETEhIlZGxP8UnUWVi4g7Wn9vy1tfPa4xW5rSBlwArAa2AY4HLoyIvYqNpAotAL4CXFp0EPXaIOBpYAYwGjgL+GlE7FJgJvXO14BdUkqjgH8EvhIR+xScSb1zAXB/0SG0UU5PKY1ofU3q6eJSlLaIGA4cA5yVUlqeUrobuAE4sdhkqkRK6ZqU0nXA34vOot5JKb2UUjo7pTQ3pbQ+pfQL4G+Af+nXiZTSrJTSqrZPW1+vKzCSeiEijgNeAH5TdBb1v1KUNmAisC6lNKfduZmAI21SFUXENuR/H2cVnUWVi4jvRcTLwKPAQuDmgiOpAhExCjgH+FTRWbTRvhYRSyLinog4sKeLy1LaRgBLO5xbClS0l5ekTRcRg4EfApenlB4tOo8ql1I6jfzn5QHANcCq7t+hGvFl4JKU0tNFB9FG+RywG7ADea22GyOi21HuspS25cCoDudGAcsKyCI1nIgYAFxJvq/09ILjaCOklNa13lqyI3Bq0XnUvYiYAhwCfLPoLNo4KaX7UkrLUkqrUkqXA/cAh3f3nrJsGD8HGBQRE1JKj7Wem4xTNFK/i4gALiE/BHR4SmlNwZG0aQbhPW314EBgF2Be/leQEcDAiNgzpfSmAnNp4yUgurugFCNtKaWXyEP650TE8IiYDhxJ/i9/1biIGBQRQ4GB5D90hkZEWf6DohFcCLweOCKltKLoMKpcRGwdEcdFxIiIGBgRhwLvA24vOpt6dDG5XE9pfV0E3AQcWmQoVSYitoiIQ9v+vouI44G3Ard2975SlLZWpwHDgMXAVcCpKSVH2urDmcAK4PPACa0fn1loIlUkInYGPkL+S2NRu/WGji84miqTyFOhzwDPA+cCn0gpXV9oKvUopfRySmlR24t8m9DKlFJL0dlUkcHkpa5agCXAx4CjUkrdrtXm3qOSJEl1oEwjbZIkSaVlaZMkSaoDljZJkqQ6YGmTJEmqA5Y2SZKkOmBpkyRJqgOWNkmSpDpgaZOkXoiIL0bEd4vOIanxWNokqXf2BB4qOoSkxmNpk6Te2Qt4sOgQkhqPpU2SuhARAyLiCxExLyIWRMRxwO7Aw0Vnk9R4BhUdQJJq2BeBtwMHAC8ANwMLU0rLCk0lqSFZ2iSpExExDvg0MDml9FTruZuAqYUGk9SwnB6VpM4dDDySUnqi3blt8CEESQWxtElS58YCi9s+iYjBwFH4EIKkgljaJKlzs4H9I2JiRIwGLgTG40ibpIJY2iSpEymlXwM/BpqB+4EWYCXwWJG5JDWuSCkVnUGSJEk9cKRNkiSpDljaJEmS6oClTZIkqQ5Y2iRJkuqApU2SJKkOWNokSZLqgKVNkiSpDljaJEmS6oClTZIkqQ78f+w7o+Na4T8jAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画出f(d)图像\n",
    "\n",
    "def f_d(d):\n",
    "    return (1/(d**2)) * (10/(3*np.sqrt(2*np.pi))) *  np.exp(-((1/d-1)**2*50)/9)\n",
    "    \n",
    "x_ref = np.arange(0.0001, 5, 0.01)\n",
    "fig = plt.figure(figsize=[10,6])\n",
    "ax = fig.add_subplot(111)\n",
    "e1, = ax.plot(x_ref, f_d(x_ref), 'b-')\n",
    "ax.set_xlabel(r'$d$', fontsize=12)\n",
    "ax.set_ylabel(r'$f(d)$', fontsize=12)\n",
    "ax.set_xlim(0, 5)\n",
    "#fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance rate: 0.099177717543844\n",
      "E(d)=1.113 \t D(d)=0.203\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAF9CAYAAACJeyWnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHa9JREFUeJzt3X+0XWdd5/H3h6a2TNNIa0NnQJsItlTTNSnDHXF0QFhV+bEGdRpdE6nUVqBQVnUQZ1WWq4Xa0sVycGRmFBjDUFt+qNiZFASdOjNaweKP4UZWcWJjsWos2ErKiiFJfzHtd/44+8rJaZK7772555z7nPdrrbNyzn7295zvXSc3/fTZz947VYUkSZLWvqdMugFJkiSdGAY7SZKkRhjsJEmSGmGwkyRJaoTBTpIkqREGO0mSpEYY7CRJkhphsJMkSWqEwU6SJKkRBjtJkqRGrJt0A5Nw1lln1ebNmyfdhiRJ0qJ27dr1YFVt7LPvTAa7zZs3Mz8/P+k2JEmSFpVkb999PRQrSZLUCIOdJElSIwx2kiRJjTDYSZIkNcJgJ0mS1AiDnSRJUiMMdpIkSY0w2EmSJDViLMEuyaGRx+NJfmFo/KIke5I8lOSOJJuGxk5JclOSLyd5IMmbRt77mLWSJEmzZCzBrqrWLzyAs4GHgVsBkpwF7ASuBc4E5oEPD5VfB5wLbAJeDFyd5KU9ayVJkmbGJA7F/gDwReD3u9cXA7ur6taqeoRBkNua5Pxu/FLghqraX1V3A+8FLutZK0mSNDMmEex+BHh/VVX3egtw18JgVR0G7gW2JDkDeMbwePd8y2K1q9a9JEnSlBprsEtyDvCdwC1Dm9cDB0Z2PQCc3o0xMr4wtljt6GdfkWQ+yfy+ffuW9wNIkiRNsXHP2F0K3FlVfzW07RCwYWS/DcDBboyR8YWxxWqPUFU7qmququY2bty4zPa1EsmTH5Ik6cSZRLC7ZWTbbmDrwoskpwHPZrB2bj9w//B493z3YrUnvHNJkqQpN7Zgl+TbgWfSnQ075DbggiTbkpwKvAX4bFXt6cbfD1yT5IzupIjXAjf3rJUkSZoZ45yx+xFgZ1UdcZi0qvYB24Abgf3A84HtQ7u8lcEJEXuBTwDvqKrbe9ZKkiTNjHz15NTZMTc3V/Pz85NuY+YcbU3dDP71kyRpSZLsqqq5PvuuW+1mpOM51gkUBj5JkpbOe8VKkiQ1wmAnSZLUCA/FalV4jTpJksbPGTtJkqRGGOwkSZIaYbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJaoTBTpIkqREGO0mSpEYY7CRJkhphsJMkSWqEwU6SJKkRBjtJkqRGGOwkSZIasW7SDUhHkzx5W9X4+5AkaS1xxk6SJKkRBjtJkqRGGOwkSZIaYbCTJElqhMFOkiSpEQY7SZKkRni5E63Y0S5NIkmSxs8ZO0mSpEYY7CRJkhphsJMkSWqEwU6SJKkRBjtJkqRGeFas1oyjnX1bNf4+JEmaVs7YSZIkNcJgJ0mS1IixBrsk25PcneRwknuTvKDbflGSPUkeSnJHkk1DNackuSnJl5M8kORNI+95zFpJkqRZMrZgl+S7gZ8FLgdOB14I/GWSs4CdwLXAmcA88OGh0uuAc4FNwIuBq5O8tHvPxWolSZJmxjhn7H4GuL6q/qiqnqiqL1TVF4CLgd1VdWtVPcIgyG1Ncn5XdylwQ1Xtr6q7gfcCl3Vji9VKkiTNjLEEuyQnAXPAxiR/keTzSX4xyVOBLcBdC/tW1WHgXmBLkjOAZwyPd8+3dM+PWbuaP48kSdI0GteM3dnAycAPAC8ALgSeC1wDrAcOjOx/gMHh2vVDr0fHWKT2CEmuSDKfZH7fvn3L/0kkSZKm1LiC3cPdn79QVfdX1YPAzwMvBw4BG0b23wAc7MYYGV8YY5HaI1TVjqqaq6q5jRs3LvsHkSRJmlZjCXZVtR/4PHC0y8nuBrYuvEhyGvBsBmvn9gP3D493z3cvVnsi+5ckSVoLxnnyxC8DP5bk6d3auTcCHwduAy5Isi3JqcBbgM9W1Z6u7v3ANUnO6E6KeC1wcze2WK0kSdLMGGewuwH4NHAPcDfwGeDGqtoHbANuBPYDzwe2D9W9lcEJEXuBTwDvqKrbAXrUSpIkzYzUDN5sc25urubn5yfdRjOOdg/XcZnBv76SpBmTZFdVzfXZ11uKSZIkNcJgJ0mS1AiDnSRJUiMMdpIkSY0w2EmSJDXCYCdJktQIg50kSVIjDHaSJEmNMNhJkiQ1wmAnSZLUCIOdJElSIwx2kiRJjTDYSZIkNcJgJ0mS1AiDnSRJUiMMdpIkSY0w2EmSJDXCYCdJktQIg50kSVIjDHaSJEmNWDfpBqSVSJ68rWr8fUiSNA2csZMkSWqEwU6SJKkRBjtJkqRGuMZOS3K0NW2SJGk6OGMnSZLUCIOdJElSIwx2kiRJjTDYSZIkNcJgJ0mS1AiDnSRJUiMMdpIkSY0w2EmSJDXCYCdJktQIg50kSVIjxhbskvxekkeSHOoefz409soke5McTvKRJGcOjZ2Z5LZubG+SV4687zFrJUmSZsm4Z+yuqqr13eM5AEm2AL8EvAo4G3gIePdQzbuAx7qxS4D3dDV9aiVJkmbGukk3wCCsfayqPgmQ5Frg7iSnA08A24ALquoQcGeS32AQ5N58vNqqOjiBn0WSJGlixj1j9/YkDyb5VJIXddu2AHct7FBV9zKYoTuvezxeVfcMvcddXc1itZIkSTNlnDN2PwX8GYPgtR34WJILgfXAgZF9DwCnA48fZ4xFao+Q5ArgCoBzzjln2T+EJEnStBrbjF1V/XFVHayqR6vqFuBTwMuBQ8CGkd03AAcXGaPH+PDn76iquaqa27hx48p+GE215MkPSZJmwSQvd1JAgN3A1oWNSZ4FnALc0z3WJTl3qG5rV8MitZIkSTNlLMEuydOSvCTJqUnWJbkEeCHw28CHgFckeUGS04DrgZ3d7N5hYCdwfZLTknwH8H3AB7q3PmbtOH4uSZKkaTKuNXYnA28Dzmewbm4P8P1V9ecASV7PIKR9HfC/gcuHat8A3AR8EfgScGVV7Qaoqt2L1EqSJM2MVNWkexi7ubm5mp+fn3Qba9JaXa82g3/NJUmNSLKrqub67OstxSRJkhphsJMkSWqEwU6SJKkRBjtJkqRGGOwkSZIaYbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJaoTBTpIkqREGO0mSpEYY7CRJkhphsJMkSWqEwU6SJKkRBjtJkqRGGOwkSZIaYbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGrFu0g1I45A8eVvV+PuQJGk1OWMnSZLUCIOdJElSIwx2kiRJjegV7JKctNqNSJIkaWX6ztjdn+Q/JZlb1W4kSZK0bH2D3cuAx4GPJbk7yU8nOWcV+5IkSdIS9Qp2VbWrqt4EPBP4CeBbgD9NckeSH01y2mo2KUmSpMUt6eSJqnoC2NM99jEIepcA9yV51YlvT5IkSX31PXnijCSvS3InsItBoLu0qs6rqouAlwD/eRX7lCRJ0iL63nni88AdDMLbR6vq0eHBqvp0ko+e6OYkSZLUX99g96yq+rvj7VBVl628HUmSJC1X3zV2lyf558MbknxrkqtXoSdJkiQtQ99g92+BPxvZ9mfAG09sO5IkSVquvsHua4CvjGx7DDh1qR+Y5NwkjyT54NC2VybZm+Rwko8kOXNo7Mwkt3Vje5O8cuT9jlkrSZI0S/oGu13AG0a2vR74k2V85ruATy+8SLIF+CXgVcDZwEPAu0f2f6wbuwR4T1fTp1aSJGlm9D154ieA/9Vdq+5e4JsYBKnvXsqHJdkO/D3wB917wCCsfayqPtntcy1wd5LTgSeAbcAFVXUIuDPJbzAIcm8+Xm1VHVxKb5IkSWtd3ztP7AbOA36OwWzbvweeU1Wj6+6OKckG4HrgJ0eGtgB3DX3WvQxm6M7rHo9X1T1D+9/V1SxWK0mSNFP6ztjRzZj96go+6wbgfVV1X5Lh7euBAyP7HgBOZ3B/2mONLVZ7hCRXAFcAnHOOt7mVJEnt6RXsknwjcCNwIYMw9Q+qatGUlORC4LuA5x5l+BCwYWTbBuAgg0OxxxpbrPYIVbUD2AEwNzdXi/UsSZK01vSdsfsVBmvrfpLBCQpL9SJgM/A33WzdeuCkJN8C3A5sXdgxybOAU4B7GAS7dUnOrarPdbtsBXZ3z3cfp1aSJGmm9A12W4DvqKonlvk5O4BfG3r97xgEvSuBpwN/mOQFDM6yvR7YuXDyQ5KdwPVJXsNgxvD7gG/v3udDx6vVyhx5xFySJE27vpc7+SRHP4zaS1U9VFUPLDwYHEJ9pKr2dSdmvJ5BSPsig/Vxw5dWeQPw1G7sV4Eruxp61EqSJM2MVC2+3CzJLwLbgZ3AA8NjVfWW1Wlt9czNzdX8/Pyk25h6rc/Y9firL0nSxCXZVVVzffbteyj2NOBjwMnANyy3MUmSJK2eXsGuqi5f7UakcTvajKSzeJKktaz3deySfDPwA8DZVXVVkucAp1TVZ1etO0mSJPXW6+SJJD/I4ASKZwKXdptPB35+lfqSJEnSEvU9K/Z64Lur6vUM7gYBg1t5bT12iSRJksapb7B7Ol+9J2sN/emKJEmSpCnRN9jtAl41sm078H9ObDuSJElarr4nT/w48D+TvBo4LclvA+cB37NqnUmSJGlJ+l7uZE+S84F/BXwcuA/4eFUdWs3mJEmS1F/vy51U1UPAr69iL5IkSVqBXsEuye9zjBMlquqFJ7QjSZIkLUvfGbv/OvL6HwOvBj54YtuRJEnScvVdY3fL6LYk/x34ZQbXuJMkSdKE9b3cydF8AfinJ6oRSZIkrUzfNXY/OrLpHwEXA390wjuSJEnSsvRdYzd6ceLDwB8A7zyx7UiSJGm5+q6xe/FqNyJJkqSV6Xso9ll99quqv1xZO5IkSVquvodi/4KvXscuHHlNu3R/FnDSCepLkiRJS9T3rNhXA78GnA+c2v35K8Crq+op3cNQJ0mSNEF9Z+xuAM6tqoe7159L8jrgHuDm1WhMkiRJS9N3xu4pwOaRbZvw0KskSdLU6Dtj907gd5P8MnAf8A3AZXi5EzUmefK2OupdkiVJmj59L3fyjiR/Cvwg8FzgfuBHq+r21WxOkiRJ/fWdsaMLcQY5SZKkKdVrjV2SU5LcmOQvkxzotn1PkqtWtz1JkiT11ffkiXcCFwCX8NVr2O0GrlyNpiRJkrR0fQ/F/mvgm6rqcJInAKrqC0meuXqtSZIkaSn6ztg9xkgITLIR+NIJ70iSJEnL0jfY3QrckuQbAZL8E+AXGdyNQpIkSVOgb7D7aeCvgT8FngZ8Dvhb4GdWpy1JkiQt1aJr7JI8BfiXwE9V1Ru7Q7APVnnZVkmSpGmy6IxdVT0BfLSqHu1e7zPUSZIkTZ++h2I/meTbVrUTSZIkrUjfYLcX+B9Jbk5yQ5LrFx59PyjJB5Pcn+TLSe5J8pqhsYuS7EnyUJI7kmwaGjslyU1d3QNJ3jTyvseslSRJmiXHDHYjd5X4WuAjDC5O/PXANww9+no7sLmqNgDfC7wtyfOSnAXsBK4FzgTmgQ8P1V0HnAtsAl4MXJ3kpV2Pi9VKkiTNjOOdPHEjg0uaALyiC2TLVlW7h192j2cDzwN2V9WtAEmuAx5Mcn5V7QEuBS6vqv3A/iTvBS5jcN/aixeplSRJmhnHC3b3JvkPDG4ddnKSy4GM7lRVN/X9sCTvZhDKngp8BvgtBgHyrqH3O5zkXmBLkr8DnjE83j3//u75lmPVAgY7SZI0U44X7LYDVwM/BJzMYOZsVAG9g11VvSHJjwH/AngR8CiwHtg3susB4PRubOH16BiL1B4hyRXAFQDnnHNO35YlSZLWjGMGu6q6B3gNQJLfqaqLTsQHVtXjwJ1Jfhi4EjgEjB7m3QAc7MYWXj8yMsYitaOfuwPYATA3N+flWiRJUnN6nRV7okLdiHUM1tjtBrYubExy2sL2bl3d/cPj3fOF9XrHrF2FfiVJkqZa38udrEiSpyfZnmR9kpOSvITBId7fBW4DLkiyLcmpwFuAzw6d/PB+4JokZyQ5H3gtcHM3tlitJEnSzBhLsGOwFu9K4PPAfuDngDdW1Uerah+wjcFJFPuB5zNY37fgrcC9DK6l9wngHVV1OwzugrFIrSRJ0szILN4dbG5urubn5yfdxtTLk86Bnk0z+CsiSZoiSXZV1VyffY93Vqwkjh5wDXuSpGk0rkOxkiRJWmUGO0mSpEYY7CRJkhphsJMkSWqEwU6SJKkRBjtJkqRGGOwkSZIaYbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJaoTBTpIkqRHrJt2AtBYlT95WNf4+JEka5oydJElSIwx2kiRJjfBQrICjH1qUJElrizN2kiRJjTDYSZIkNcJgJ0mS1AiDnSRJUiMMdpIkSY0w2EmSJDXCYCdJktQIg50kSVIjDHaSJEmNMNhJkiQ1wmAnSZLUCIOdJElSI9ZNugGpFcmTt1WNvw9J0uxyxk6SJKkRBjtJkqRGGOwkSZIaYbCTJElqxFiCXZJTkrwvyd4kB5N8JsnLhsYvSrInyUNJ7kiyaaT2piRfTvJAkjeNvPcxayVJkmbJuGbs1gH3Ad8JfC1wLfDrSTYnOQvY2W07E5gHPjxUex1wLrAJeDFwdZKXAvSolSRJmhljudxJVR1mENAWfDzJXwHPA74O2F1VtwIkuQ54MMn5VbUHuBS4vKr2A/uTvBe4DLgduHiRWkmSpJkxkTV2Sc4GzgN2A1uAuxbGuhB4L7AlyRnAM4bHu+dbuufHrF3N/iVJkqbR2INdkpOBDwG3dLNq64EDI7sdAE7vxhgZXxhjkdrRz70iyXyS+X379q3sh5AkSZpCYw12SZ4CfAB4DLiq23wI2DCy6wbgYDfGyPjC2GK1R6iqHVU1V1VzGzduXPbPIEmSNK3GFuySBHgfcDawraq+0g3tBrYO7Xca8GwGa+f2A/cPj3fPdy9Wu0o/hiRJ0tQa54zde4BvBl5RVQ8Pbb8NuCDJtiSnAm8BPjt08sP7gWuSnJHkfOC1wM09ayVJkmbGuK5jtwl4HXAh8ECSQ93jkqraB2wDbgT2A88Htg+Vv5XBCRF7gU8A76iq2wF61EoTlRz9IUnSakhVTbqHsZubm6v5+flJtzFVDBvjNYO/dpKkZUqyq6rm+uzrLcUkSZIaYbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJasS6STcgzaKjXV7GS6BIklbKGTtJkqRGGOwkSZIaYbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJaoTBTpIkqREGO0mSpEYY7CRJkhrhLcWkKeFtxiRJK+WMnSRJUiMMdpIkSY0w2EmSJDXCYCdJktQIg50kSVIjDHaSJEmN8HIn0hTzEiiSpKVwxk6SJKkRBjtJkqRGeCh2Bh3t8J4kSVr7nLGTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJaoRnxUprjBctliQdy9hm7JJclWQ+yaNJbh4ZuyjJniQPJbkjyaahsVOS3JTky0keSPKmvrWSJEmzZJyHYv8WeBtw0/DGJGcBO4FrgTOBeeDDQ7tcB5wLbAJeDFyd5KU9ayVJkmbG2IJdVe2sqo8AXxoZuhjYXVW3VtUjDILc1iTnd+OXAjdU1f6quht4L3BZz1pJkqSZMQ0nT2wB7lp4UVWHgXuBLUnOAJ4xPN4937JY7Sr3LEmSNHWmIditBw6MbDsAnN6NMTK+MLZY7RGSXNGt8Zvft2/fipuWJEmaNtMQ7A4BG0a2bQAOdmOMjC+MLVZ7hKraUVVzVTW3cePGFTctTZPkyQ9J0uyZhmC3G9i68CLJacCzGayd2w/cPzzePd+9WO0q9yxJkjR1xnm5k3VJTgVOAk5KcmqSdcBtwAVJtnXjbwE+W1V7utL3A9ckOaM7KeK1wM3d2GK1kiRJM2OcM3bXAA8DbwZ+uHt+TVXtA7YBNwL7gecD24fq3srghIi9wCeAd1TV7QA9aiVJkmZGagYvWT83N1fz8/OTbmNiXH81G2bwV1uSmpRkV1XN9dl3GtbYSZIk6QQw2EmSJDVi3aQbkLQ6jnbI3cOzktQ2Z+wkSZIaYbCTJElqhMFOkiSpEa6xk2aI6+4kqW3O2EmSJDXCYCdJktQIg50kSVIjDHaSJEmN8OQJacb1vXewJ1lI0vRzxk6SJKkRBjtJkqRGGOwkSZIa4Ro7Sb0cay2ea+8kaXo4YydJktQIg50kSVIjPBTbuL6XspCWy/vPStL0MNhJOuEMe5I0GR6KlSRJaoTBTpIkqREGO0mSpEa4xk7SWLjuTpJWnzN2kiRJjXDGTtLE9L0cjzN7ktSPM3aSJEmNcMZO0tRzZk+S+nHGTpIkqRHO2ElqhmfeSpp1BjtJTTPsSZolBjtJM6fvmr1jMRhKmlausZMkSWqEM3aStEQrnfEb5QygpBPFYNeQE/0fG0nj4TpASSdKE4dik5yZ5LYkh5PsTfLKSfckSSuRPPmxkv1W+jmS1oZWZuzeBTwGnA1cCPxmkruqavdk25KkE6dv6BpXOHOmUZo+az7YJTkN2AZcUFWHgDuT/AbwKuDNE21OktaglQTDtTrjZyBVK9Z8sAPOAx6vqnuGtt0FfOeE+hmLtfqPpyRNI/9N1WKWEv4nOZvdQrBbDxwY2XYAOH14Q5IrgCu6l48m+b9j6E2r4yzgwUk3oWXz+1u7/O7WNr+/FVhp+F9h/XP67thCsDsEbBjZtgE4OLyhqnYAOwCSzFfV3Hja04nm97e2+f2tXX53a5vf39qVZL7vvi2cFXsPsC7JuUPbtgKeOCFJkmbKmg92VXUY2Alcn+S0JN8BfB/wgcl2JkmSNF5rPth13gA8Ffgi8KvAlYtc6mTHWLrSavH7W9v8/tYuv7u1ze9v7er93aU8x1uSJKkJrczYSZIkzTyDnSRJUiNmKth5T9m1K8lVSeaTPJrk5kn3o6VJckqS93W/dweTfCbJyybdl/pJ8sEk9yf5cpJ7krxm0j1p6ZKcm+SRJB+cdC/qJ8nvdd/Zoe7x54vVzFSw48h7yl4CvCfJlsm2pJ7+FngbcNOkG9GyrAPuY3BHmK8FrgV+PcnmCfak/t4ObK6qDcD3Am9L8rwJ96Slexfw6Uk3oSW7qqrWd49FL1Q8M8Fu6J6y11bVoaq6E1i4p6ymXFXtrKqPAF+adC9auqo6XFXXVdVfV9UTVfVx4K8Aw8EaUFW7q+rRhZfd49kTbElLlGQ78PfA70y6F62umQl2HPuess7YSWOW5GwGv5NeSHyNSPLuJA8Be4D7gd+acEvqKckG4HrgJyfdi5bl7UkeTPKpJC9abOdZCna97ikraXUlORn4EHBLVe2ZdD/qp6rewODfyxcwuCj8o8ev0BS5AXhfVd036Ua0ZD8FPAt4JoNr2X0syXFny2cp2PW6p6yk1ZPkKQzuCvMYcNWE29ESVdXj3TKWrweunHQ/WlySC4HvAt456V60dFX1x1V1sKoerapbgE8BLz9ezbrxtDYV/uGeslX1uW6b95SVxiRJgPcxOHnp5VX1lQm3pOVbh2vs1ooXAZuBvxn8CrIeOCnJt1TVP5tgX1qeAnK8HWZmxs57yq5tSdYlORU4icE/SqcmmaX/MWnBe4BvBl5RVQ9Puhn1k+TpSbYnWZ/kpCQvAX4I+N1J96ZedjAI4Rd2j/8C/Cbwkkk2pcUleVqSlyz89y7JJcALgd8+Xt3MBLvOUu8pq+lxDfAw8Gbgh7vn10y0I/WWZBPwOgb/YXlg6JpMl0y4NS2uGBx2/TywH/g54I1V9dGJdqVequqhqnpg4cFgWdIjVbVv0r1pUSczuMzXPuBB4MeA76+q417LznvFSpIkNWLWZuwkSZKaZbCTJElqhMFOkiSpEQY7SZKkRhjsJEmSGmGwkyRJaoTBTpJOsCQ3J3nbpPuQNHsMdpIkSY0w2EmSJDXCYCdJK5TkuUn+JMnBJB8GTp10T5Jmk8FOklYgydcAHwE+AJwJ3Apsm2hTkmaWwU6SVubbGNys+z9W1Veq6r8Bn55wT5JmlMFOklbmGcAXqqqGtu2dVDOSZpvBTpJW5n7gmUkytO2cSTUjabYZ7CRpZf4Q+H/AjydZl+Ri4Fsn3JOkGWWwk6QVqKrHgIuBy4D9wL8Bdk6yJ0mzK0cuC5EkSdJa5YydJElSIwx2kiRJjTDYSZIkNcJgJ0mS1AiDnSRJUiMMdpIkSY0w2EmSJDXCYCdJktQIg50kSVIj/j98ibo2NE7jlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 使用拒绝抽样得到d的分布\n",
    "\n",
    "N = 100000\n",
    "i = 0\n",
    "k = 0\n",
    "x = []\n",
    "while i<N:\n",
    "    x0 = np.random.uniform(0, 5)\n",
    "    y0 = np.random.rand()*2\n",
    "    y1 = f_d(x0)\n",
    "    if y0 <= y1:\n",
    "        x.append(x0)\n",
    "        i += 1\n",
    "    k += 1\n",
    "\n",
    "fig = plt.figure(figsize=[10,6])\n",
    "ax = fig.add_subplot(111)\n",
    "ax.hist(x, bins=100, color='b')\n",
    "ax.set_xlabel(r'd', fontsize=12)\n",
    "ax.set_ylabel(r'frequency', fontsize=12)\n",
    "ax.set_xlim(0, 5)\n",
    "#fig.show()\n",
    "\n",
    "print ('Acceptance rate:', N/np.float(k))\n",
    "print('E(d)=%(m).3f \\t D(d)=%(d).3f' %{'m':np.mean(x), 'd':np.var(x)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\n\\n采样范围不同时得到不同结果\\n0-5:   1.112, 0.206\\n0-10:  1.124, 0.279\\n0-50:  1.140, 0.538\\n0-100: 1.141, 0.728\\n0-150: 1.149, 1.322\\n0-200: 1.151, 1.941\\n\\n'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "\n",
    "采样范围不同时得到不同结果\n",
    "0-5:   1.112, 0.206\n",
    "0-10:  1.124, 0.279\n",
    "0-50:  1.140, 0.538\n",
    "0-100: 1.141, 0.728\n",
    "0-150: 1.149, 1.322\n",
    "0-200: 1.151, 1.941\n",
    "\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_0(m)=1.109\n",
      "D_0(d)=0.206\n"
     ]
    }
   ],
   "source": [
    "# 使用scipy计算\n",
    "\n",
    "def f_d_e(d):\n",
    "    return (1/d) * (10/(3*np.sqrt(2*np.pi))) *  np.exp(-((1/d-1)**2*50)/9)\n",
    "\n",
    "E = integrate.quad(f_d_e, 0, 5)[0]\n",
    "print ('E_0(m)=%(m).3f' %{'m':E})\n",
    "\n",
    "def f_d_d(d):\n",
    "    return ((d - E)**2) * (1/(d**2)) * (10/(3*np.sqrt(2*np.pi))) *  np.exp(-((1/d-1)**2*50)/9)\n",
    "\n",
    "D = integrate.quad(f_d_d, 0, 5)[0]\n",
    "print ('D_0(d)=%(d).3f' %{'d':D})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(theta)=0.999 \t D(theta)=0.090\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAF+CAYAAAALGoc0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XuQXGl55/nvU1VSS6VLS2pdunWpqkwB27h7aS7aIIIdcHvxeGkHY7M0QbTp5dqDPCbweI1nvUQ0t8WWd2wvG7trD8Ztq6fpbuEBG/CFHTNrZsCsicFhcWk8CkM7U6pS6Ybu15JUkurdPyqPSFVXlaoqM+vkyfx+IjJUmee8mU9eVPrpPed9MlJKSJIkqbh68i5AkiRJjTHQSZIkFZyBTpIkqeAMdJIkSQVnoJMkSSo4A50kSVLBGegkSZIKbk6BLiLeHxF7I+JqRDw1ZdvrI+L7ETEWEV+NiMG6bXdExJMRcT4ijkXEBxZjrCRJUjeZ6wzdEeDXgSfrb4yI9cAXgA8D64C9wGfrdvkY8GJgEPgJ4Fcj4g2LMFaSJKlrxHy+KSIifh3YmlJ6V+36TuBdKaXX1K6vAE4Cr0gpfT8iDgPvTin9v7Xtvwa8OKX0SCvHNvqiSJIkFUmj59DdBzyXXUkpXQKqwH0RsRbYXL+99vN9rRw7tcCI2Fk7XLy3FgQlSZI6Sl+D41cCJ6bcdg5YVduWXZ+6rZVjb5FSegJ4AmD9+vVpx44dvz/z05EkSWoP3/rWt06mlDbMZd9GA91FYPWU21YDF2rbsutXpmxr5dgZDQ0NsXfv3tl2kSRJagsRMTLXfRs95LoPeKDugVcA24F9KaUzwNH67bWf97VybIPPR5IkqXDm2rakLyKWAb1Ab0Qsi4g+4IvA/RHxcG37R4Dv1S1MeBr4UESsjYh7gfcCT9W2tXKsJElS15jrDN2HgMvAB4H/sfbzh1JKJ4CHgV3AGeDVwCN14z7K5GKFEeCvgd9OKX0ZoMVjJUmSusa82pYU3Y4dO5Ln0EmSpCKIiG+llHbMZV+/+kuSJKngDHSSJEkFZ6CTJEkqOAOdJElSwRnoJEmSCs5AJ0ltYs+ePQwNDdHT08PQ0BB79uzJuyRJBdHoV39Jkppgz5497Ny5k7GxMQBGRkbYuXMnAI8++miepUkqAGfoJKkNPP744zfDXGZsbIzHH388p4okFYmBTpLawMGDB+d1uyTVM9BJUhsYGBiY1+2SVM9AJ0ltYNeuXSxbtuyW2/r7+9m1a1dOFUkqEgOdJLWBRx99lMcee+zm9YGBAZ544gkXREiaE1e5SlKbuPvuu2/+/PWvf53BwcEcq5FUJM7QSVKbqFarN3/ev39/jpVIKhoDnSS1iUqlwrZt2wA4cOBAztVIKhIDnSS1iWq1yoMPPkhvb68zdJLmxUAnSW3g0qVLHD16lHvvvZeBgQEDnaR5MdBJUhvIAtz27dspl8secpU0LwY6SWoDlUoFgBe96EWUy2Vn6CTNi4FOktpAtsJ1+/btlEoljh8/zsWLF3OuSlJRGOgkqQ1UKhXuuusu1qxZQ7lcBlzpKmnuDHSS1Aaq1Srbt28HMNBJmjcDnSS1gUqlwote9CLgR4HO8+gkzZWBTpJyNj4+zsGDB2/O0K1bt45Vq1YZ6CTNmYFOknI2PDzMxMTEzRm6iLB1iaR5MdBJUs7qV7hmbF0iaT4MdJKUs/oedJlSqcT+/ftJKeVVlqQCMdBJUs6q1SorVqxg48aNN28rl8tcuXKFY8eO5ViZpKIw0ElSzrIVrhFx8zZbl0iaDwOdJOWsvgddxtYlkubDQCdJObpx4wb79++/5fw5gMHBQcBAJ2luDHSSlKPDhw8zPj7+ghm6ZcuWsWXLFg+5SpoTA50k5Wi6Fa4ZW5dImisDnSTlaLoedJmsdYkk3Y6BTpJyVKlUWLp0KVu3bn3BtnK5zOHDh7l69WoOlUkqEgOdJOWoWq1SKpXo7e19wbZyuUxKiZGRkRwqk1QkBjpJylHWg246pVIJcKWrpNsz0ElSTlJK0/agy9iLTtJcGegkKSfHjx/n4sWLM87Q3X333SxbtszWJZJuy0AnSTmZbYUrQE9PjytdJc2JgU6ScjJbD7qMgU7SXBjoJCkn1WqVnp4ehoaGZtwnay6cUlq8wiQVjoFOknJSqVQYGBhg6dKlM+5TLpc5f/48Z86cWcTKJBWNgU6ScjLbCteMrUskzYWBTpJyMlsPuoytSyTNhYFOknJw9uxZTp06NecZOluXSJqNgU6ScpC1LLndDN2qVatYv369M3SSZmWgk6Qc3K4HXb1spaskzcRAJ0k5yHrQzTXQechV0mwMdJKUg2q1yt13382KFStuu2+5XGZkZITr168vQmWSishAJ0k5mMsK10ypVOL69escOnSoxVVJKioDnSTlYC496DJZ6xIPu0qaiYFOkhbZ5cuXOXz48Jxn6OxFJ+l2Gg50ETEUEf8+Is5ExLGI+N2I6Ktte3lEfCsixmp/vrxuXETEb0bEqdrltyIi6rYveKwktbMsmM11hm7r1q309vYa6CTNqBkzdJ8EjgP3AC8Hfhx4X0QsBf4MeBZYC3wa+LPa7QA7gTcBDwAvA94I/DxAI2Mlqd1lK1znOkPX19fH4OCggU7SjJoR6ErA51JKV1JKx4AvA/cBDwJ9wP+ZUrqaUvq/gQD+u9q4dwKfSCkdSikdBj4BvKu2rZGxktTW5tODLmPrEkmzaUag+7+ARyKiPyK2AA/xo1D3vZRSqtv3e7Xbqf35XN2256ZsW+jYW0TEzojYGxF7T5w4Me8nJ0nNVqlUWLNmDevWrZvzGJsLS5pNMwLdXzMZps4Dh4C9wJ8CK4FzU/Y9B6yq/Tx1+zlgZe1cuEbG3iKl9ERKaUdKaceGDRvm+dQkqfmq1eqcD7dmSqUSJ06c4OLFiy2qSlKRNRToIqIH+A/AF4AVwHomz3n7TeAisHrKkNXAhdrPU7evBi7WZuUaGStJba1SqczrcCvYukTS7BqdoVsHbAN+t3au2yng3wI/DewDXjZl1uxltdup/flA3bYHpmxb6FhJalvXrl1jZGRk3jN0ti6RNJuGAl1K6SRwAPiFiOiLiDVMLlh4DvgacAP4lxFxR0S8vzbsP9X+fBr4QERsiYjNwK8AT9W2NTJWktrWyMgIN27cmPcMXalUAgx0kqbXjHPo3gy8ATgBVIDrwC+nlMaZbC3yDuAs8B7gTbXbAX4f+Avg74H/Avw/tdtoZKwktbNshet8Z+jWrVvH6tWrPeQqaVp9jd5BSum7TLYZmW7bd4BXzbAtAb9auzR1rCS1q6wH3Xxn6CLCla6SZuRXf0nSIqpWqyxfvpx77rln3mNLpZKBTtK0DHSStIiq1Srbt29nId9WmDUXdkG/pKkMdJK0iCqVyrzPn8uUy2WuXLnCsWPHmlyVpKIz0EnSIpmYmGD//v3zPn8uY+sSSTMx0EnSIjly5AhXrlxZ8AydrUskzcRAJ0mLJGtZstAZusHBQSLC1iWSXsBAJ0mLJGtZstAZumXLlrFlyxZn6CS9gIFOkhZJtVqlr6+Pbdu2Lfg+bF0iaToGOklaJJVKhVKpRF/fwnu621xY0nQMdJK0SLIedI0ol8s3F1dIUsZAJ0mLIKXUUA+6TKlUIqXEyMhIkyqT1AkMdJK0CE6dOsX58+ebMkMHti6RdCsDnSQtgkZXuGayQGfrEkn1DHSStAga7UGXufvuu1m2bJkzdJJuYaCTpEVQqVSIiJvf9rBQ2X0Y6CTVM9BJ0iKoVqts3bqVZcuWNXxf5XLZQ66SbmGgk6RF0IwVrpmsF11KqSn3J6n4DHSStAia0YMuUyqVOH/+PKdPn27K/UkqPgOdJLXYhQsXOH78eFNn6MDWJZJ+xEAnSS3WrBWuGVuXSJrKQCdJLdasHnSZbKWsM3SSMgY6SWqxZs/QrVy5kg0bNhjoJN1koJOkFqtUKmzcuJFVq1Y17T5tXSKpnoFOklqsmStcM1nrEkkCA50ktVwze9BlSqUSIyMjXL9+van3K6mYDHSS1EJXrlzh0KFDLZmhu3HjBqOjo029X0nFZKCTpBY6cOAAKaWmz9DZukRSPQOdJLVQs1e4ZmxdIqmegU6SWqjZPegyW7dupa+vz0AnCTDQSVJLVatVVq9ezV133dXU++3r62NwcNBDrpIAA50ktVS2wjUimn7fti6RlDHQSVILtaIHXaZUKhnoJAEGOklqmevXr3PgwIGmnz+XKZfLnDx5kgsXLrTk/iUVh4FOklpkdHSU69evt2yGztYlkjIGOklqkVatcM3YukRSxkAnSS3Sqh50mWyGzkAnyUAnSS1SqVRYtmwZmzdvbsn9r127ljvvvNNDrpIMdJLUKtVqlXK5TE9Pa37VRoQrXSUBBjpJapmsB10r2YtOEhjoJKklUkot7UGXKZfLDA8PMzEx0dLHkdTeDHSS1AJHjx7l8uXLizJDd+XKFY4dO9bSx5HU3gx0ktQCrV7hmrF1iSQw0ElSS7S6B13G1iWSwEAnSS1RrVbp7e1lYGCgpY8zODhIRNi6ROpyBjpJaoFKpcLQ0BBLlixp6ePccccdbNmyxRk6qcsZ6CSpBRZjhWvG1iWSDHSS1AKL0YMuUy6XPeQqdTkDnSQ12enTpzl79uyizdCVSiUOHz7MlStXFuXxJLUfA50kNdlirXDNZCtdh4eHF+XxJLUfA50kNdli9aDLZIHOw65S9zLQSVKTZTN0WdBqNXvRSTLQSVKTVatVtmzZwvLlyxfl8TZt2sTy5csNdFIXa1qgi4hHIuIfIuJSRFQj4rW1218fEd+PiLGI+GpEDNaNuSMinoyI8xFxLCI+MOU+FzxWkvKymCtcASKCUqlkoJO6WFMCXUT8U+A3gXcDq4DXAfsjYj3wBeDDwDpgL/DZuqEfA14MDAI/AfxqRLyhdp8LHitJeVrMHnQZW5dI3a1ZM3T/K/DxlNI3U0oTKaXDKaXDwJuBfSmlP04pXWEyhD0QEffWxr0D+LWU0pmU0j8AfwC8q7atkbGSlIuLFy9y7NixRZ2hA27O0KWUFvVxJbWHhgNdRPQCO4ANEVGJiEMR8bsRsRy4D3gu2zeldAmoAvdFxFpgc/322s/31X5uZGx9fTsjYm9E7D1x4kSjT1eSZpUd9sxjhu7ChQucOnVqUR9XUntoxgzdJmAJ8BbgtcDLgVcAHwJWAuem7H+OycOyK+uuT91Gg2NvSik9kVLakVLasWHDhrk/K0lagMXuQZexdYnU3ZoR6C7X/vydlNLRlNJJ4P8Afhq4CKyesv9q4EJtG1O2Z9tocKwk5WKxe9BlbF0idbeGA11K6QxwCJjuxI19wAPZlYhYAWxn8ty4M8DR+u21n/c1Yawk5aJSqbB+/XruvPPORX3coaEhwEAndatmLYr4t8AvRsTG2vlt/xPwJeCLwP0R8XBELAM+AnwvpfT92ringQ9FxNraYof3Ak/VtjUyVpJykccKV4CVK1eyceNGD7lKXapZge7XgL8Dngf+AfgOsCuldAJ4GNgFnAFeDTxSN+6jTC50GAH+GvjtlNKXARoZK0l5WewedPXK5bIzdFKX6mvGnaSUrgHvq12mbvsKcO8LBk1uuwq8p3aZbvuCx0rSYrt69Sqjo6O5zNDBZOuSb37zm7k8tqR8+dVfktQkw8PDTExM5DpDd/DgQa5fv57L40vKj4FOkpokrxWumXK5zI0bNxgdHc3l8SXlx0AnSU2SVw+6TKlUAlzpKnUjA50kNUm1WmXlypXk1cTcXnRS9zLQSVKTZCtcIyKXx9+6dSt9fX22LpG6kIFOkpokrx50md7eXoaGhpyhk7qQgU6SmuDGjRvs378/t/PnMqVSyUAndSEDnSQ1waFDh7h27VquM3Rgc2GpWxnoJKkJ8l7hmimXy5w6dYrz58/nWoekxWWgk6QmyLsHXSZrXeLCCKm7GOgkqQkqlQp33HEHW7duzbUOW5dI3clAJ0lNUK1WKZVK9PTk+2s1C3TO0EndxUAnSU2Q9aDL29q1a7nzzjudoZO6jIFOkhqUUsq9B109V7pK3cdAJ0kN+uEPf8ilS5faYoYOJgOdh1yl7mKgk6QGtcsK10wW6CYmJvIuRdIiMdBJUoPapQddplQqcfXqVY4ePZp3KZIWiYFOkhpUrVbp6elhcHAw71IAW5dI3chAJ0kNqlQqDA4OsnTp0rxLAWxdInUjA50kNaidVrgCDAwMEBHO0EldxEAnSQ1qlx50mewbKwx0Uvcw0ElSA86cOcPp06fbaoYObF0idRsDnSQ1IGtZ0k4zdGBzYanbGOgkqQHt1oMuUyqVOHLkCJcvX867FEmLwEAnSQ3IetBlK0vbRVbP8PBwvoVIWhQGOklqQLVa5Z577mHFihV5l3ILW5dI3cVAJ0kNaLcVrplSqQTYXFjqFgY6SWpAu/Wgy2zatInly5cb6KQuYaCTpAUaGxvjyJEjbTlDFxG2LpG6iIFOkhYom/1qxxk6mDzs6gyd1B0MdJK0QNkK13acoYMf9aJLKeVdiqQWM9BJ0gK1aw+6TLlc5uLFi5w6dSrvUiS1mIFOkhaoUqmwbt061q5dm3cp08pal3jYVep8BjpJWqB2XeGasXWJ1D0MdJK0QO3agy5joJO6h4FOkhZgfHyckZGRtp6hW7FiBZs2bbJ1idQFDHSStAAjIyNMTEy09Qwd2LpE6hYGOklagHZf4ZrJWpdI6mwGOklagHbvQZcpl8uMjo5y7dq1vEuR1EIGOklagGq1evMctXZWKpW4ceMGo6OjeZciqYUMdJK0AJVKhe3btxMReZcyK3vRSd3BQCdJC9DuPegyBjqpOxjoJGmeJiYm2L9/f9ufPwewZcsWlixZYusSqcMZ6CRpng4fPszVq1cLMUPX29vL4OCgM3RShzPQSdI8FWWFa8bWJVLnM9BJ0jwVpQddplwue8hV6nAGOkmap0qlwpIlS9i2bVvepcxJqVTi1KlTnDt3Lu9SJLWIgU6S5qlarVIqlejt7c27lDnJVro6Syd1LgOdJM1TpVIpzPlzYKCTuoGBTpLmIaVUmB50GXvRSZ3PQCdJ83DixAkuXLhQqBm6NWvWsGbNGgOd1MEMdJI0D0Vb4ZqxdYnU2ZoW6CLixRFxJSKerbvtbRExEhGXIuJPI2Jd3bZ1EfHF2raRiHjblPtb8FhJapWi9aDL2LpE6mzNnKH7N8DfZVci4j7g94G3A5uAMeCTU/Yfr217FPi92piGxkpSK1WrVSKCoaGhvEuZl1KpxIEDB5iYmMi7FEkt0JRAFxGPAGeB/1h386PAX6SUvp5Sugh8GHhzRKyKiBXAw8CHU0oXU0p/A/w5kwGu0bGS1DKVSoWBgQHuuOOOvEuZl3K5zPj4OEeOHMm7FEkt0HCgi4jVwMeBX5my6T7guexKSqnK5KzaS2qXGyml5+v2f642ptGxU+vbGRF7I2LviRMn5v8EJalO0Va4ZmxdInW2ZszQ/RqwO6U0OuX2lcDUtuTngFW32dbo2FuklJ5IKe1IKe3YsGHDbZ6KJM2uaD3oMqVSCbB1idSp+hoZHBEvB34SeMU0my8Cq6fcthq4AEzMsq3RsZLUEufOnePkyZOFnKEbHBwkIgx0UodqKNABDwJDwMGIgMnZs96I+DHgy8AD2Y4RUQbuAJ5nMpT1RcSLU0r/WNvlAWBf7ed9DYyVpJbIWpYUcYZu6dKlbNu2zUAndahGA90TwL+ru/6vmAx4vwBsBP5zRLwW+DaT59l9IaV0ASAivgB8PCL+OfBy4GeB19TuZ08DYyWpJYragy5j6xKpczV0Dl1KaSyldCy7MHmo9EpK6URKaR/wL5gMZ8eZPMftfXXD3wcsr237I+AXamNoZKwktUrWg66oga5UKjlDJ3WoRmfobpFS+tiU658BPjPDvqeBN81yXwseK0mtUK1W2bRpEytXrsy7lAUpl8scPXqUy5cvs3z58rzLkdREfvWXJM1RUVe4ZrLWJcPDw/kWIqnpDHSSNEdF7UGXsXWJ1LkMdJI0B5cvX+bQoUMdMUNnoJM6j4FOkuYgWx1a5Bm6jRs30t/f70pXqQMZ6CRpDrIVrkWeoYsIV7pKHcpAJ0lzUPQedJlyuWygkzqQgU6S5qBSqbBmzRrWrVuXdykNyQJdSinvUiQ1kYFOkuYgW+Fa+5rDwiqXy1y6dImTJ0/mXYqkJjLQSdIcFL0HXcbWJVJnMtBJ0m1cu3aNkZGRwp8/B7YukTqVgU6SbuPgwYNcv369o2bobF0idRYDnSTdRqescAXo7+9n06ZNztBJHcZAJ0m30Qk96OrZukTqPAY6SbqNarXK8uXLueeee/IupSnK5bKHXKUOY6CTpNuoVCqUy+XCtyzJlEolDh48yLVr1/IuRVKTGOgk6Taq1WrHHG6FyRm6iYkJDh48mHcpkprEQCdJs5iYmLjZVLhT2LpE6jwGOkmaxdGjR7ly5UrHzdCBrUukTmKgk6RZZCtcO2mGbvPmzSxZssQZOqmDGOgkaRZZD7pOmqHr7e1laGjIQCd1EAOdJM2iUqnQ19fHwMBA3qU0la1LpM5ioJOkWVSrVYaGhujr68u7lKYqlUrO0EkdxEAnSbOoVCoddf5cplwuc/r0ac6ePZt3KZKawEAnSTNIKXVcD7qMK12lzmKgk6QZnDp1inPnznXsDB0Y6KROYaCTpBl04grXTKlUAmwuLHUKA50kzaATe9Bl1qxZw9q1aw10Uocw0EnSDKrVKhFx8/Bkp7F1idQ5DHSSNINKpcKWLVtYtmxZ3qW0hK1LpM5hoJOkGXTqCtdMuVxmeHiYiYmJvEuR1CADnSTNoFN70GXK5TLj4+McOXIk71IkNchAJ0nTuHDhAsePH+/oGTpXukqdw0AnSdPIWpZ0+gwdGOikTmCgk6RpdHIPuszAwAA9PT0GOqkDGOgkaRqd3IMus3TpUrZt22brEqkDGOgkaRrVapUNGzawevXqvEtpKVuXSJ3BQCdJ0+j0Fa6ZcrlsoJM6gIFOkqbR6T3oMuVymWPHjjE2NpZ3KZIaYKCTpCmuXr3K6OhoV8zQZa1LhoeH8y1EUkMMdJI0xYEDB0gpdc0MHdi6RCo6A50kTdENK1wzBjqpMxjoJGmKbuhBl9mwYQP9/f22LpEKzkAnSVNUKhVWrVrF+vXr8y6l5SLCla5SBzDQSdIU2QrXiMi7lEVhoJOKz0AnSVN0Sw+6TLlcvrkQRFIxGegkqc7169cZHh7uivPnMqVSiUuXLnHixIm8S5G0QAY6SaozOjrKtWvXum6GDlzpKhWZgU6S6nTTCtdMFuhc6SoVl4FOkup0Uw+6zNDQEOAMnVRkBjpJqlOtVrnjjjvYsmVL3qUsmv7+fu6++24DnVRgBjpJqlOpVCiXy/T0dNevR1uXSMXWXb+xJOk2sh503SZrXSKpmAx0klSTUqJarXbV+XOZUqnE6Ogo4+PjeZciaQEaDnQRcUdE7I6IkYi4EBHfiYiH6ra/PiK+HxFjEfHViBicMvbJiDgfEcci4gNT7nvBYyVpvo4dO8bY2FjXztBNTExw8ODBvEuRtADNmKHrA0aBHwfuBD4MfC4ihiJiPfCF2m3rgL3AZ+vGfgx4MTAI/ATwqxHxBoBGxkrSQnTjCteMrUukYutr9A5SSpeYDFeZL0XEAeBVwF3AvpTSHwNExMeAkxFxb0rp+8A7gHenlM4AZyLiD4B3AV8G3tzAWEmat27sQZcplUqArUukomr6OXQRsQl4CbAPuA94LttWC39V4L6IWAtsrt9e+/m+2s+NjK2vZ2dE7I2IvX6tjaTZVCoVent7GRwcvP3OHWbz5s0sXbrUQCcVVFMDXUQsAfYAn67Noq0Ezk3Z7RywqraNKduzbTQ49qaU0hMppR0ppR0bNmyY3xOS1FWq1SqDg4MsWbIk71IWXW9vL0NDQwY6qaCaFugiogd4BhgH3l+7+SKwesquq4ELtW1M2Z5ta3SsJM1bpVLpyvPnMqVSyXPopIJqSqCLiAB2A5uAh1NK12qb9gEP1O23AtjO5LlxZ4Cj9dtrP+9rwlhJmrdu7UGXsbmwVFzNmqH7PeClwD9LKV2uu/2LwP0R8XBELAM+AnyvdjgW4GngQxGxNiLuBd4LPNWEsZI0L6dPn+bMmTNdPUNXLpc5c+YMZ8+ezbsUSfPUjD50g8DPAy8HjkXExdrl0ZTSCeBhYBdwBng18Ejd8I8yudBhBPhr4LdTSl8GaGSsJM1XN69wzdi6RCquZrQtGQFilu1fAe6dYdtV4D21S1PHStJ8dHMPukx965JXvOIVOVcjaT786i9J4kczdNksVTfKnrvn0UnFY6CTJCZn6DZv3kx/f3/epeTmzjvvZN26dR5ylQrIQCdJuMI1UyqVnKGTCshAJ0nYgy5j6xKpmAx0krrepUuXOHbsmDN0TAa64eFhbty4kXcpkubBQCep62ULIpyhmzzkeu3aNY4cOZJ3KZLmwUAnqevZg+5HXOkqFZOBTlLXswfdjxjopGIy0EnqetVqlbvuuos1a9bkXUruBgYG6OnpsXWJVDAGOkldzxWuP7JkyRK2bdvmDJ1UMAY6SV3PHnS3snWJVDwGOkldbXx8nIMHDzpDV6dcLnvIVSoYA52krjY8PMzExIQzdHVKpRLHjh1jbGws71IkzZGBTlJXc4XrC2UrXZ2lk4rDQCepq9mD7oVsXSIVj4FOUlerVCqsWLGCjRs35l1K23CGTioeA52krpatcI2IvEtpG+vXr2fFihXO0EkFYqCT1NXsQfdCEWHrEqlgDHSSutaNGzc4cOCA589Nw9YlUrEY6CR1rUOHDjE+Pu4M3TRKpRL79+8npZR3KZLmwEAnqWu5wnVm5XKZsbExjh8/nncpkubAQCepa9mDbma2LpGKxUAnqWtVq1WWLl3K1q1b8y6l7ZRKJcDWJVJRGOgkda1KpUKpVKK3tzfvUtrO0NAQ4AydVBQGOkldK+tBpxfq7+/nnnvuMdBJBWGgk9SVUkr2oLsNW5dIxWGgk9SVjh8/zqVLl5yhm0XWukRS+zPQSepKrnC9vXK5zOjoKOPj43mXIuk2DHSSupI96G6vXC6TUuLgwYN5lyLpNgx0krpSpVIl6LQ+AAANbElEQVShp6fn5mpOvVDWusTDrlL7M9BJ6krVapWBgQGWLl2adylty+bCUnEY6CR1JVe43t7mzZtZunSpgU4qAAOdpK5kD7rbyw5J27pEan8GOkld5+zZs5w6dcoZujkol8vO0EkFYKCT1HVc4Tp3BjqpGAx0krqOPejmrlwuc/bsWc6cOZN3KZJmYaCT1HWyGToD3e1lrUs8j05qbwY6SV1lz549/MZv/AYA9913H3v27Mm5ovZm6xKpGPryLkCSFsuePXvYuXMnY2NjAIyMjLBz504AHn300TxLa1s2F5aKwRk6SV3jgx/84M0wlxkbG+Pxxx/PqaL2d+edd7Ju3ToPuUptzkAnqaNNTEzwla98hZ/7uZ/j0KFD0+7jd5XOzpWuUvvzkKukjjQ6OspTTz3Fk08+yfDwMGvWrGHlypVcvHjxBfsODAzkUGFxlMtlvv3tb+ddhqRZOEMnqWOMj4/zJ3/yJzz00EMMDg7ykY98hHK5zJ49ezhy5Aif+tSn6O/vv2VMf38/u3btyqniYiiXy4yMjHDjxo28S5E0A2foJBXevn372L17N8888wwnT55ky5YtPP7447z73e++uUoTfrTw4fHHH+fgwYMMDAywa9cuF0TcRqlU4tq1axw+fNjZTKlNGegkFdL58+f57Gc/y+7du/nbv/1blixZws/8zM/w2GOP8VM/9VP09vZOO+7RRx81wM1TfesSA53Ungx0kgojpcQ3vvENdu/ezec+9znGxsb4sR/7MT7xiU/w9re/nQ0bNuRdYkfKAt2BAwd48MEH8y1G0rQMdJLa3rFjx3j66ad58skn+cEPfsDKlSt529vexmOPPcarX/1qIiLvEjvatm3b6OnpcaWr1MYMdJLa0vXr1/nLv/xLdu/ezZe+9CVu3LjBa17zGnbv3s1b3/pWVq5cmXeJXWPJkiUMDAwY6KQ2ZqCT1Fb+8R//kSeffJJPf/rTHD16lI0bN/LLv/zLvOc97+GlL31p3uV1LXvRSe3NQCcpd2NjY3z+85/nD//wD/n6179OT08PDz30EI899hhvfOMbWbJkSd4ldr1SqcSXvvSlvMuQNAMDnaRcpJT41re+xe7du/nMZz7D+fPn2b59O7t27eKd73wnW7ZsybtE1SmXy/zwhz/k0qVLrFixIu9yJE1hoJO0qE6dOsWePXvYvXs33/ve91i2bBlvectbeOyxx3jd615HT4/9zttR/UrX+++/P+dqJE1V6N+cEbEuIr4YEZciYiQi3pZnPXv27GFoaIienh6GhobYs2dPnuXcwtoWxtoWZmptzz77LH/1V3/FI488wubNm/mlX/ollixZwic/+UmOHj3KM888w4MPPmiYa2PPP/88AC972cva/vNmbXNjbQvTtrWllAp7Af4I+CywEvgnwDngvpn2f9WrXpVa5dlnn039/f0JuHnp7+9Pzz77bMse09qsrSi1RUQC0tq1a9Mv/uIvpu9+97t5l6l5ePbZZ9Py5csL83mzNmvrlNqAvWmOmSgm9y+eiFgBnAHuTyk9X7vtGeBwSumD043ZsWNH2rt3b0vqGRoaYmRk5AW3L126lFe+8pUtecy5+va3v834+PgLbre22VnbwsxU2/r16xkdHWXZsmU5VKVGzPT7rbe3l82bNxMRN3sBZj834/pc9nnuuecK93fB2mZXxNoGBwcZHh5u+uNFxLdSSjvmsm+Rz6F7CXAjC3M1zwE/Xr9TROwEdgIt/cqagwcPTnv7+Pg4q1evbtnjzsV0H77sdmubmbUtzEy1nTp1yjBXUDP9frtx4wavf/3rpx45mfH6bNsWer2IfxesbXZFrG2mvyOLaq5Tee12AV4LHJty23uBr800ppWHXAcHB2+Zgs0ug4ODLXtMa7M2a9NiaOf31NqsrZNrYx6HXIt8BvJFYGpUXw1cyKEWdu3aRX9//y239ff3s2vXrjzKuYW1LYy1LUw716aFaef31NoWxtoWpp1ry32mbaEXYAUwDry47rangX8905hWztClNHmy5ODgYIqINDg42BYncGasbWGsbWHauTYtTDu/p9a2MNa2MItZG92wKAIgIv4dk9Od/xx4OfDvgdeklPZNt38rF0VIkiQ103wWRRT5kCvA+4DlwHEmW5j8wkxhTpIkqVMVeZUrKaXTwJvyrkOSJClPRZ+hkyRJ6noGOkmSpIIz0EmSJBWcgU6SJKngDHSSJEkFZ6CTJEkqOAOdJElSwRnoJEmSCq7QX/01XxFxAhiZZtN64OQil9NufA18DcDXAHwNwNcAfA3A1wDyfw0GU0ob5rJjVwW6mUTE3rl+V1qn8jXwNQBfA/A1AF8D8DUAXwMo1mvgIVdJkqSCM9BJkiQVnIFu0hN5F9AGfA18DcDXAHwNwNcAfA3A1wAK9Bp4Dp0kSVLBOUMnSZJUcAY6SZKkgjPQSZIkFVxXBLqIeH9E7I2IqxHx1DTbXx8R34+IsYj4akQMznJfQ7V9xmpjfrKlxbdIRFyccrkREb8zw77vqm2v3//BRS656SLiaxFxpe45/WCWfSMifjMiTtUuvxURsZj1NltE3BERuyNiJCIuRMR3IuKhWfbvmM9BRKyLiC9GxKXa83/bDPt13PsO83vvO+l9n2quvwM6+HPQlf8OzJYJipwHuiLQAUeAXweenLohItYDXwA+DKwD9gKfneW+/gj4DnAX8DjwJxExpy7O7SSltDK7AJuAy8AfzzLkP9ePSSl9bVEKbb331z2n/2qW/XYCbwIeAF4GvBH4+cUosIX6gFHgx4E7mfw78LmIGJplTKd8Dv4NMM7kZ/9R4Pci4r5p9uvE9x3m/953yvs+nbn8DujIz0EX/zswbSYoeh7oikCXUvpCSulPgVPTbH4zsC+l9McppSvAx4AHIuLeqTtGxEuAVwIfTSldTil9Hvh74OHWVb8o3gIcB/6/vAtpY+8EPpFSOpRSOgx8AnhXviU1JqV0KaX0sZTScEppIqX0JeAA8Kq8a2uliFjB5N/ZD6eULqaU/gb4c+Dt0+zece87dO9734CO/BxM0TX/DsySCQqdB7oi0N3GfcBz2ZWU0iWgWrt9un33p5Qu1N323Az7Fsk7gafT7D1sXhERJyPi+Yj4cET0LVZxLfa/1Z7XN25z+OCWzwmd8b7fIiI2AS8B9s2yWyd8Dl4C3EgpPV9320zvZ8e/7zCn974T3veZzOV3QDd8Drr534FMofOAgQ5WAuem3HYOWNXgvoUQEQNMHnb59Cy7fR24H9jI5P8+fg74n1tfXcv9L0AZ2MJk88i/iIjtM+w79b0/B6zshPNoACJiCbAH+HRK6fsz7NYpn4NG/s531PsOc3rvO+V9n85cfwd09Oegy/8dqFfoPFD4QFc7qTXNcPmbOdzFRWD1lNtWAxca3Dc383xN3gH8TUrpwEz3l1Lan1I6UDs08/fAx5mcnm9bc3kNUkp/m1K6kFK6mlL6NPAN4KdnuMup7/1q4OJt/jebq7l+DiKiB3iGyXPK3j/T/RXxczCDRv7Ot/37Ph9zee876H1/gXn8DujozwEd+u/AAhQ6DxQ+0KWUHkwpxQyXfzKHu9jH5ImuwM3za7Yz/aGHfUA5IuoT+AMz7Jubeb4m72D2/5VN+xBAW//PdIGfi9me1y2fE9rwfZ9qLq9BbYZhN5MnRD+cUro2n4egzT8HM3ge6IuIF9fdNtP7Wbj3fa4aeO+L+r7PxUzPrWM/BzUd+e/AAhQ6DxQ+0M1FRPRFxDKgF+iNiGV1x/6/CNwfEQ/X9vkI8L3pDj3Uzrn5LvDR2n38D0yuePr84jyT5oqI1zB5qGG2VU1ExEO1c2yonRz6YeDPWl9h60TEmoj477PPQkQ8CrwO+A8zDHka+EBEbImIzcCvAE8tUrmt9HvAS4F/llK6PNuOnfI5qJ0X8wXg4xGxIiL+W+BnmZypmqpT33eY43vfKe/7VPP8HdCxn4Nu/HdglkxQ7DyQUur4C5MrVdKUy8fqtv8k8H0ml2x/DRiq2/Yp4FN114dq+1wGfgD8ZN7Pr4HX5feBZ6a5fYDJ6eSB2vX/HfghcAnYz+RU+5K862/wuW8A/o7J6fGzwDeBf1q3/bVMHlLJrgfwW8Dp2uW3qH0XclEvwGDt78KV2vudXR7t9M8Bky0J/rT2XA4Cb+uW9/12730nv+9TXoMZfwd0y+eg9ty67t8BZskEFDgPRK0oSZIkFVRXHHKVJEnqZAY6SZKkgjPQSZIkFZyBTpIkqeAMdJIkSQVnoJMkSSo4A50kSVLBGegkSZIKzkAnSZJUcAY6SVqAiNgeEacj4pW165sj4mREPJhzaZK6kF/9JUkLFBHvBT4AvIrJL/b++5TSv8q3KkndyEAnSQ2IiD8HSkx+wfd/k1K6mnNJkrqQh1wlqTF/ANwP/I5hTlJenKGTpAWKiJXAc8BXgYeA/zqldDrfqiR1IwOdJC1QROwGVqWU3hoRTwBrUkpvzbsuSd3HQ66StAAR8bPAG4B/UbvpA8ArI+LR/KqS1K2coZMkSSo4Z+gkSZIKzkAnSZJUcAY6SZKkgjPQSZIkFZyBTpIkqeAMdJIkSQVnoJMkSSo4A50kSVLB/f8U9cv+ANNDiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 使用随机抽样求theta的分布\n",
    "\n",
    "mu = 1.0\n",
    "sigma2 = 0.09\n",
    "sigma = np.sqrt(sigma2)\n",
    "x1 = np.random.normal(loc=mu, scale=sigma, size=100000)\n",
    "\n",
    "xgrid = np.arange(-10,12,1.5)\n",
    "xcenter = (xgrid[1:]+xgrid[:-1])/2\n",
    "hx, xedge = np.histogram(x1, xgrid)\n",
    "\n",
    "fig = plt.figure(figsize=[10,6])\n",
    "ax = fig.add_subplot(111)\n",
    "e1, = ax.plot(xcenter, hx, 'ko-')\n",
    "ax.set_xlabel('x', fontsize=12)\n",
    "#fig.show()\n",
    "\n",
    "print('E(theta)=%(m).3f \\t D(theta)=%(d).3f' %{'m':np.mean(x1), 'd':np.var(x1)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_1(d)=1.166 \t D_1(d)=8.913\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAF9CAYAAACJeyWnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGOtJREFUeJzt3X+Q7WddH/D3J7kx0SQXSXNJJWIiNoiGMaFu1alFsagIo9UmOo2kUH5oECZagQ46DiCCjGOtohalxvJLUSu0gQhaaEdRi9WWjU7AK2kQNQZM6g0TL7kJJAif/nHOxc02e/e7e+/u2X3O6zVzZs95nu+z+zlz7t77vs/3+zzf6u4AALD/nbboAgAAODUEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBHFh0AYtw/vnn98UXX7zoMgAANnXjjTfe2d2Hphy7lMHu4osvzurq6qLLAADYVFXdOvVYp2IBAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIMQ7AAABiHYAQAMQrADABiEYAcAMAjBDgBgEAcWXQBLpmrjvu7dqwMABmTGDgBgEIIdAMAgnIplZ5zolCsAsCPM2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIKyKZe/YaCWtjYsBYBIzdgAAgxDsAAAGIdgBAAxCsAMAGIRgBwAwCMEOAGAQgh0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBC7Fuyq6uKq+o2ququq7qiqV1XVgXnf5VV1Y1XdO/96+ZpxVVU/WlUfmT/+bVXVmv4NxwIALJPdnLH72SR/neRzklye5KuTPLeqPiPJDUnemOShSd6Q5IZ5e5Jck+RbklyW5EuSfGOSZyfJhLEAAEtjN4Pd5yd5U3d/vLvvSPKOJJcmeXySA0l+srvv6+6fTlJJ/ul83L9K8uPd/aHu/nCSH0/y9HnfZmMBAJbGbga7n0pyVVV9VlVdmORJ+btw997u7jXHvnfenvnXm9b03bSu70RjAQCWxm4Gu9/JLHB9NMmHkqwmeWuSc5IcXXfs0STnzp+v7z+a5Jz5dXabjf20qrqmqlaravXIkSMn+VYAAPaeXQl2VXVakncmuT7J2UnOz+yauB9NcizJwXVDDia5e/58ff/BJMfms3Sbjf207r6uu1e6e+XQoUMn94YAAPag3ZqxOy/JI5K8an4t3EeSvC7Jk5McTvIla1e6ZrZI4vD8+eHMFk4cd9m6vhONBQBYGrsS7Lr7ziR/nuQ5VXWgqj47s0URNyX57SSfTPI9VXVmVV07H/Zb86+/kOT5VXVhVT08yQuSvH7et9lYAIClsZvX2F2R5BuSHEnyp0n+Nsnzuvv+zLYzeVqSv0nyzCTfMm9Pkp9L8rYk70vyx0l+fd6WCWMBAJZGPXBB6XJYWVnp1dXVRZcxtgecHT9JS/hnFACOq6obu3tlyrFuKQYAMIgDiy6Afe5UzswBACfFjB0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBuKUYe9+JblvWvXt1AMAeZ8YOAGAQgh0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIMQ7AAABiHYAQAMQrADABiEYAcAMAjBDgBgEIIdAMAgBDsAgEEIdgAAgxDsAAAGIdgBAAxCsAMAGMSBRRcAJ6Xqwdu7d7cOANgDzNgBAAxCsAMAGIRgBwAwCMEOAGAQgh0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIMQ7AAABiHYAQAMQrADABiEYAcAMAjBDgBgEIIdAMAgBDsAgEEIdgAAgxDsAAAGsavBrqquqqr3V9U9VfXBqnrcvP0JVXVzVd1bVe+qqovWjDmzql5bVR+tqjuq6vnrvueGYwEAlsmuBbuq+rokP5rkGUnOTfJVSf6sqs5Pcn2SFyc5L8lqkl9dM/SlSS5JclGSr0nywqr6hvn33GwsAMDS2M0Zux9K8rLu/oPu/lR3f7i7P5zkiiSHu/vN3f3xzILcZVX16Pm4pyV5eXff1d3vT/LzSZ4+79tsLADA0tiVYFdVpydZSXKoqv60qj5UVa+qqs9McmmSm44f2933JPlgkkur6qFJHr62f/780vnzDcfu5PsBANiLdmvG7oIkZyT51iSPS3J5kscmeVGSc5IcXXf80cxO156z5vX6vmwy9gGq6pqqWq2q1SNHjmz/nQAA7FG7Few+Nv/677v79u6+M8lPJHlykmNJDq47/mCSu+d9Wdd/vC+bjH2A7r6uu1e6e+XQoUPbfiMAAHvVrgS77r4ryYeS9IN0H05y2fEXVXV2ki/I7Nq5u5LcvrZ//vzwZmNPZf0AAPvBbi6eeF2S766qh82vnfveJG9P8pYkj6mqK6vqrCQvSfLe7r55Pu4Xkryoqh46XxTxnUleP+/bbCwAwNLYzWD38iTvSXJLkvcn+aMkr+juI0muTPKKJHcl+fIkV60Z94OZLYi4NcnvJPmx7n5HkkwYy7Kq2vgBAIOq7gc7Ozq2lZWVXl1dXXQZ+8doYWgJ/8wDsH9V1Y3dvTLlWLcUAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIMQ7AAABiHYAQAMQrADABiEYAcAMAjBDgBgEIIdAMAgBDsAgEEIdgAAg5gU7Krq9J0uBACAkzN1xu72qvqpqlrZ0WoAANi2qcHuSUk+meRtVfX+qvqBqvq8HawLAIAtmhTsuvvG7n5+kguTPC/JFyd5X1W9q6qeWVVn72SRAABsbkuLJ7r7U0lunj+OZBb0rk5yW1U99dSXBwDAVFMXTzy0qp5dVe9OcmNmge5p3f2o7n5Ckicm+ekdrBMAgE0cmHjch5K8K7PwdkN337e2s7vfU1U3nOriAACYbmqwe2R3/98THdDdTz/5cgAA2K6p19g9o6r+0dqGqvqyqnrhDtQEAMA2TA12/zrJn6xr+5Mk33tqywEAYLumBrvPSPKJdW33Jznr1JYDAMB2TQ12NyZ57rq270ryh6e2HAAAtmvq4onnJfnv873qPpjkHyS5IMnX7VRhAABszaRg192Hq+pRSb4pyecmuT7J27v72E4WBzui6sHbu3e3DgA4xabO2GUe4n5lB2sBAOAkTAp2VfX5SV6R5PIk56zt6+7P24G6AADYoqkzdr+c2bV1L0hy786VAwDAdk0Ndpcm+cru/tROFgMAwPZN3e7kd5M8dicLAQDg5EydsfuLJO+squuT3LG2o7tfcqqLAgBg66YGu7OTvC3JGUkesXPlAACwXVP3sXvGThcCAMDJmbyPXVV9UZJvTXJBd19bVV+Y5Mzufu+OVQcAwGSTFk9U1bdltoDiwiRPmzefm+QndqguAAC2aOqq2Jcl+bru/q4kn5y33ZTksh2pCgCALZsa7B6WWZBLkl7z1c01AQD2iKnB7sYkT13XdlWS/31qywEAYLumLp74niT/raqeleTsqnpnkkcl+fodqwwAgC2Zut3JzVX16CTfmOTtSW5L8vbuPraTxQEAMN3k7U66+94kb9rBWgAAOAmTgl1V/Y9ssFCiu7/qlFYEAMC2TJ2x+4/rXv/9JM9K8sZTWw4AANs19Rq7N6xvq6r/kuR1me1xBwDAgk3d7uTBfDjJl5yqQgAAODlTr7F75rqmz0pyRZI/OOUVAQCwLVOvsVu/OfE9Sf5nklee2nIAANiuqdfYfc1OFwIAwMmZeir2kVOO6+4/O7lyAADYrqmnYv80f7ePXeWBe9rV/GsnOf0U1QUAwBZNXRX7rCT/Kcmjk5w1//rLSZ7V3afNH0Id+1vVxg8A2Aemzti9PMkl3f2x+esPVNWzk9yS5PU7URgAAFszdcbutCQXr2u7KE69AgDsGVNn7F6Z5Leq6nVJbkvyiCRPj+1OAAD2jKnbnfxYVb0vybcleWyS25M8s7vfsZPFAQAw3dQZu8xDnCAHALBHTbrGrqrOrKpXVNWfVdXRedvXV9W1O1seAABTTV088cokj0lydf5uD7vDSZ6zE0UBALB1U4PdP0/ylO7+/SSfSpLu/nCSC7f6A6vqkqr6eFW9cU3bU6rq1qq6p6reWlXnrek7r6reMu+7taqesu77bTgWAGCZTA1292fd9XhVdSjJR7bxM38myXvWfJ9Lk/xckqcmuSDJvUl+dt3x98/7rk7y6vmYKWMBAJbG1GD35iRvqKrPT5Kq+pwkr8rsbhSTVdVVSf4myW+uab46ydu6+3e7+1iSFye5oqrOraqzk1yZ5MXdfay7353k1zILciccu5W6AABGMDXY/UCSv0jyviSfneQDSf4qyQ9N/UFVdTDJy5K8YF3XpUluOv6iuz+Y2Qzdo+aPT3b3LWuOv2k+ZrOxAABLZdNgV1WnJfknSb6vu8/J7JTnud39vO6+fws/6+VJXtPdt61rPyfJ0XVtR5Ocu0nfZmPXv49rqmq1qlaPHDmyhbIBAPaHTYNdd38qyQ3dfd/89ZHu7k2GPUBVXZ7ka/Pgd6o4luTguraDSe7epG+zsQ/Q3dd190p3rxw6dGgr5QMA7AtTNyj+3ar6iu7+g23+nMdndq/Zv6yqZDbTdnpVfXFmmx5fdvzAqnpkkjOT3JLZCtwDVXVJd39gfshlmW21kvnXjcYCACyVqcHu1iT/tapuyOxesZ+esevul0wYf10euNDi32QW9J6T5GFJfr+qHpfkDzO7Du/67r47Sarq+iQvq6rvSHJ5km9O8o/n3+eXTjQWAGCZbHgqdt1dJR6S5K2ZBbrPTfKINY9Ndfe93X3H8Udmp1A/Pj+tezjJd2UW0v46s+vjnrtm+HOTfOa871eSPGc+JhPGAgAsjdrocrmqOtrdD5k//2h3r7+Wbd9aWVnp1dXVRZexf8xOny+3rV1WCgCnTFXd2N0rU4490anYD1bVj2d2HdsZVfWMJP/fv/Dd/drtlQkAwKl0omB3VZIXJvn2JGckedqDHNNJBDsAgD1gw2A33xT4O5Kkqn6zu5+wa1UBALBlk+48IdQBAOx9U28pBgDAHjd1HzuWgdWvALCvmbEDABiEYAcAMAjBDgBgEK6xgyk2uv7QHSkA2EPM2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIMQ7AAABiHYAQAMQrADABiEYAcAMAjBDgBgEIIdAMAgBDsAgEEIdgAAgxDsAAAGIdgBAAziwKILgH2t6sHbu3e3DgCIGTsAgGEIdgAAgxDsAAAGIdgBAAxCsAMAGIRgBwAwCMEOAGAQgh0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIMQ7AAABnFg0QXAkKo27uvevToAWCpm7AAABiHYAQAMQrADABiEYAcAMAjBDgBgEIIdAMAgBDsAgEEIdgAAgxDsAAAGIdgBAAxCsAMAGIR7xcJu2+g+su4hC8BJ2pUZu6o6s6peU1W3VtXdVfVHVfWkNf1PqKqbq+reqnpXVV20buxrq+qjVXVHVT1/3ffecCwAwDLZrVOxB5LcluSrkzwkyYuTvKmqLq6q85NcP287L8lqkl9dM/alSS5JclGSr0nywqr6hiSZMBYAYGnsyqnY7r4ns4B23Nur6s+TfGmSv5fkcHe/OUmq6qVJ7qyqR3f3zUmeluQZ3X1Xkruq6ueTPD3JO5JcsclYAIClsZDFE1V1QZJHJTmc5NIkNx3vm4fADya5tKoemuTha/vnzy+dP99w7E7WDwCwF+16sKuqM5L8UpI3zGfVzklydN1hR5OcO+/Luv7jfdlk7Pqfe01VrVbV6pEjR07uTQAA7EG7Guyq6rQkv5jk/iTXzpuPJTm47tCDSe6e92Vd//G+zcY+QHdf190r3b1y6NChbb8HAIC9ateCXVVVktckuSDJld39iXnX4SSXrTnu7CRfkNm1c3cluX1t//z54c3G7tDbAADYs3Zzxu7VSb4oyTd198fWtL8lyWOq6sqqOivJS5K8d83ih19I8qKqemhVPTrJdyZ5/cSxAABLY7f2sbsoybOTXJ7kjqo6Nn9c3d1HklyZ5BVJ7kry5UmuWjP8BzNbEHFrkt9J8mPd/Y4kmTAWAGBpVC/hbvcrKyu9urq66DL2no3uiMDuWMLfRQA2V1U3dvfKlGPdKxYAYBCCHQDAIAQ7AIBB7MotxYAJTnSNo+vvAJjAjB0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBA2KF42J9oEl71ro8/NxsUArGHGDgBgEIIdAMAgBDsAgEEIdgAAgxDsAAAGYVUs7GcnWuVsxSzA0jFjBwAwCMEOAGAQgh0AwCAEOwCAQQh2AACDEOwAAAYh2AEADMI+djCqjfa4s78dwLDM2AEADEKwAwAYhGAHADAIwQ4AYBCCHQDAIAQ7AIBBCHYAAIOwjx0sG/vbAQzLjB0AwCAEOwCAQQh2AACDcI0dMLPRtXeJ6+8A9gkzdgAAgxDsAAAG4VQssDlbpADsC2bsAAAGIdgBAAzCqVhg+6ykBdhTzNgBAAxCsAMAGIRTscDOsJIWYNeZsQMAGIQZO2B3WXABsGMEu1Gd6B9PAGBIgh2wd7guD+CkuMYOAGAQZuyAvW87lxaY5QOWkBk7AIBBCHYAAINwKhYYk9O3wBIS7ACOsyoX2OcEO4DN2FQZ2CcEO4CT4ZQvsIcIdgC7zSlfYIcMsSq2qs6rqrdU1T1VdWtVPWXRNQFsWdXWH9v5XsCwRpmx+5kk9ye5IMnlSX69qm7q7sOLLQtgh20nqO1WuDMDCbtu3we7qjo7yZVJHtPdx5K8u6p+LclTk3z/Qovbaf7nDexlu/F31F4Ij06ts4fs+2CX5FFJPtndt6xpuynJVy+oHgB2y17+D+5erm0kJwrQe/Uz2MHQP0KwOyfJ0XVtR5Ocu7ahqq5Jcs385X1V9ce7UBs74/wkdy66CLbN57d/+ez2tzE/v70a3k5k6zV/4dQDRwh2x5IcXNd2MMndaxu6+7ok1yVJVa1298rulMep5vPb33x++5fPbn/z+e1fVbU69dgRVsXekuRAVV2ypu2yJBZOAABLZd8Hu+6+J8n1SV5WVWdX1Vcm+eYkv7jYygAAdte+D3Zzz03ymUn+OsmvJHnOJludXLcrVbFTfH77m89v//LZ7W8+v/1r8mdXbTk2AMAQRpmxAwBYeoIdAMAglirYuafs/lZV11bValXdV1WvX3Q9TFdVZ1bVa+a/d3dX1R9V1ZMWXRfTVNUbq+r2qvpoVd1SVd+x6JrYuqq6pKo+XlVvXHQtTFdVvz3/3I7NH//nRMcvVbDLA+8pe3WSV1fVpYstiS34qyQ/nOS1iy6ELTuQ5LbM7gjzkCQvTvKmqrp4gTUx3Y8kubi7Dyb5Z0l+uKq+dME1sXU/k+Q9iy6Cbbm2u8+ZP064WfHSBLs195R9cXcf6+53Jzl+T1n2ge6+vrvfmuQji66Frenue7r7pd39F939qe5+e5I/TyIc7APdfbi77zv+cv74ggWWxBZV1VVJ/ibJby66FnbW0gS7bHxPWTN2sMuq6oLMfidtJL5PVNXPVtW9SW5OcnuS31hwSUxUVQeTvCzJCxZdC9v2I1V1Z1X9XlU9/kQHLlOwm3RPWWBnVdUZSX4pyRu6++ZF18M03f3czP6+fFxmm8Lfd+IR7CEvT/Ka7r5t0YWwLd+X5JFJLsxsP7u3VdWGM+bLFOwm3VMW2DlVdVpmd4W5P8m1Cy6HLeruT84vY/ncJM9ZdD1srqouT/K1SV656FrYnu7+X919d3ff191vSPJ7SZ680fEHdq+0hfv0PWW7+wPzNveUhV1SVZXkNZktXnpyd39iwSWxfQfiGrv94vFJLk7yl7NfwZyT5PSq+uLu/ocLrIvt6yS1UefSzNi5p+z+V1UHquqsJKdn9hfTWVW1TP852e9eneSLknxTd39s0cUwTVU9rKquqqpzqur0qnpikm9P8luLro1JrssshF8+f/yHJL+e5ImLLIppquqzq+qJx/+9q6qrk3xVknduNGZpgt3cVu8py97yoiQfS/L9Sf7l/PmLFloRk1TVRUmendk/LHes2Y/p6gWXxuY6s9OuH0pyV5J/l+R7u/uGhVbFJN19b3ffcfyR2WVJH+/uI4uujUnOyGybryNJ7kzy3Um+pbs33MvOvWIBAAaxbDN2AADDEuwAAAYh2AEADEKwAwAYhGAHADAIwQ4AYBCCHcApVlWvr6ofXnQdwPIR7AAABiHYAQAMQrADOElV9diq+sOquruqfjXJWYuuCVhOgh3ASaiqz0jy1iS/mOS8JG9OcuVCiwKWlmAHcHK+IrMbdf9kd3+iu/9zkvcsuCZgSQl2ACfn4Uk+3N29pu3WRRUDLDfBDuDk3J7kwqqqNW2ft6higOUm2AGcnN9P8rdJvqeqDlTVFUm+bME1AUtKsAM4Cd19f5Irkjw9yV1J/kWS6xdZE7C86oGXhQAAsF+ZsQMAGIRgBwAwCMEOAGAQgh0AwCAEOwCAQQh2AACDEOwAAAYh2AEADEKwAwAYxP8DfNjN2AcVipIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 尝试直接使用theta抽样值得到d的分布\n",
    "\n",
    "x1 = x1[x1>0]\n",
    "d1 = 1/x1\n",
    "\n",
    "fig = plt.figure(figsize=[10,6])\n",
    "ax = fig.add_subplot(111)\n",
    "ax.hist(d1, bins=10000, color='r')\n",
    "ax.set_xlabel(r'd', fontsize=12)\n",
    "ax.set_ylabel(r'frequency', fontsize=12)\n",
    "ax.set_xlim(0, 5)\n",
    "#fig.show()\n",
    "\n",
    "print('E_1(d)=%(m).3f \\t D_1(d)=%(d).3f' %{'m':np.mean(d1), 'd':np.var(d1)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
